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SUMMARY 

This  document  represents  the  final  report  under  Contract  No.  F49620-83-K- 
0032-P00001  from  the  Air  Force  Office  of  Scientific  Research  to  Virginia 
Polytechnic  Institute  and  State  University.  To  accommodate  the  fact  that  the 
Principal  Investigator  (J.L.  Junkins)  accepted  a  position  at  Texas  A&M 
University,  effective  September  1,  1985,  the  effort  was  performed  at  Texas  A&M 
under  Sub-Contract  No.  4174-352142-1  from  Virginia  Polytechnic  Institute. 

Significant  progress  is  reported  on  methodology  to  optimize  open  and 
closed  loop  control  laws  for  flexible  vehicles.  Also  a  new  method  for 
simutaneous  structure/controller  design  optimization  is  reported. 

DISCUSSION  OF  RESULTS 

Since  the  effort  in  this  contract  resulted  in  ten  manuscripts  which 
document  the  results  in  detail,  we  append  here  to  these  manuscripts  and 
provide  below  a  brief  statement  of  the  main  results  contained  in  these  papers. 
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Nonlinear  Feedback  Control 

In  Attachment  1,  we  present  a  generalization  of  classical  linear- 
quadratic  regulator  theory  to  accownodate  polynomial  nonlinearities  in 
dynamical  models  of  the  form 


Xi=  a^Xj  +  cijkxjx)<  + - +  b^Uj  (i=l,2,...,n)  (1) 

where  x^  are  state  variables,  u^  are  control  variables,  and  summation  over 
repeated  indices  is  implied.  We  consider  polynomial  nonlinear  feedback  of  the 


form 


For  the  case  of  a  quadratic  performance  index,  we  show  in  Attachment  No.  1 
that  the  gains  satisfy  a  sequence  of  differential  equations.  The  linear  gains 
(d^j)  are  determined  by  solving  the  usual  Riccati  equation,  whereas  the 
quadratic  gains  (e1jk),  cubic  gains  (f i )  and  all  higher  order  gains  satisfy 
linear  differential  equations.  Numerical  implementations  have  been  carried 
out  and  the  validity  of  the  formulation  has  been  established.  Several 
applications  are  described  in  Attachment  1  which  shows  that  the  nonlinear 
feedback  terms  are  significant  and  constructive  in  designing  optimal  nonlinear 
controls.  In  all  cases  studied,  except  one,  convergence  has  been  reliably 
achieved.  This  convergence  failure  was  found  to  depend  upon  the  weights 
selected  in  the  performance  index  and  was  easily  eliminated.  More  generally, 
conditions  which  guarantee  convergence  remain  a  difficult  unresolved  issue 
which  requires  further  research.  The  most  significant  practical  difficulty 
associated  with  this  approach  lies  in  the  system  specific  algebra  required  to 
derive  the  differential  equations  satisfied  by  the  higher  order  gains.  We  are 
investigating  the  use  of  algebraic  manipulators  (MACSYMA  and  SMP)  to  automate 
the  derivation  and  coding  of  these  equations. 

In  Attachment  2,  we  extend  these  ideas  by  using  canonical  state  variables 
(the  euler  parameters  and  the  corresponding  conjugate  momentum  variables). 
This  formulation  is  very  elegant  for  spacecraft  attitude  maneuvers.  While  we 
clearly  establish  the  validity  of  the  formulations  (through  analytical  and 
numerical  comparisons  to  the  results  in  Attachment  1),  the  advantages  of  this 
canonical  variable  approach  have  not  been  established.  Also,  we  consider  in 
Attachment  2  a  Liapunov  approach  to  design  of  nonlinear  feedback  control; 
these  results  look  very  promising,  especially  for  systems  of  moderate 
dimensions. 


Nonlinear  Open  Loop  Control 
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In  Attachments  3,4,5,  we  present  a  novel  method  for  computing  nonlinear 
open  loop  spacecraft  maneuvers  by  a  perturbation  method.  As  with  the  problems 
addressed  in  Attachments  1  and  2,  we  consider  only  polynomial 
nonlinearities.  We  establish  a  quasi-analytical ,  non-iterative  method  (based 
upon  asymptotic  expansions)  which  we  demonstrate  to  work  well  on  several 
problems.  The  method  presented  is  attractive  because  it  appears  suitable  for 
semi -automation  and  also  because  it  has  worked  well  on  most  examples  tried  to 
date. 


Structural  Identification 

In  Attachment  6,  we  present  a  time-domain  method  for  structural 
identification  which  combines  both  free  and  forced  response  measurements  to 
determine  system  mass  and  stiffness  estimates.  The  method  is  shown  to  work 
well  for  "academic"  structures  of  moderate  dimensionality.  However 
difficulties  associated  with  lack  of  uniqueness,  rank  deficiency  and  data 
requirements  are  noted  which  limit  the  practical  utility  of  this  approach. 
Recent  research  has  established  a  new  frequency  domain  approach  which  promises 
to  circumvent  many  (if  not  most)  of  these  difficulties.  The  method  is  based 
upon  parameterization  of  the  frequency  response  function  in  terms  of  physical 
system  parameters  and  recovering  estimates  of  these  parameters  to  fit  the 
frequency  response  in  a  least  square  sense.  Rank  deficiencies  are  addressed 
using  a  singular  value  decomposition  algorithm.  The  details  of  this  method 
will  be  presented  in  N.G.  Creamer's  forthcoming  dissertation,  including 
applications  to  structures  with  linear  viscoelastic  models. 
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Structure/Controller  Optimization 

In  Attachments  7  and  8,  we  present  homotopy  methods  for  simultaneous 
optimization  of  structural  design  parameters,  sensor/actuator  locations,  and 
the  feedback  gains  in  an  output  feedback  control  law.  The  algorithm  of 
Attachment  7  is  based  upon  a  "minimum  modification"  strategy  which  is  shown 
convergent  with  as  many  as  60  design  variables.  The  algorithm  of  Attachment  8 
is  based  upon  sequential  linear  programing.  This  approach  is  especially 
well-suited  to  high  dimensioned  problems  involving  a  large  number  of 
inequality  constraints.  In  Attachment  8,  we  consider  a  simple  example  with 
eigenvalue  placement  constraints  and  two  alternative  optimization  criteria 
(minimum  mass,  minimum  closed-loop  eigenvalue  sensitivity);  both  structure  and 
control  parameters  are  simultaneously  iterated.  Several  variations  in  problem 
statement  and  starting  iterative  support  the  validity  and  usefulness  of  both 
the  minimum  correction  and  sequential  linear  programming  homotopy  algorithms. 

In  Attachments  9  and  10,  we  extend  the  methodology  of  Attachments  7  and  8 
to  consider  tuning  of  optimal  quadratic  regulators.  In  Attachment  10,  we 
establish  that  the  weight  matrices  selected  for  the  usual  quadratic 
performance  index  have  a  strong  Impact  upon  eigenvalue  placement  and  other 
stability/performance  robustness  measures.  We  also  establish  an  algorithm 
which  we  have  found  useful  for  "optimal  tuning"  of  the  weight  matrices;  we 
have  successfully  iterated  as  many  as  150  weight  elements  to  optimize 
eigenvalue  placement  and  robustness  indices. 

The  most  significant  unifying  feature  of  the  methods  in  attachments  7-10 
Is  the  use  of  homotopy  methods.  This  approach  is  most  Important  In  practical 
situations  for  which  the  initially  stated  contraints  have  no  feasible 
solution.  The  homotopy  method  gradually  imposes  constraints,  and  as  a 
consequence,  convergence  failures  are  informative.  Conflicting  and  active 


constraints  can  be  easily  identified,  leading  to  "least  compromised"  revisions 
of  the  problem  statement.  This  feature  appears  to  be  a  key  ingredient  in 
developing  practical  optimization  strategies  for  high-dimensioned  systems. 
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Optimal  Nonlinear  Feedack  Control 
for  Spacecraft  Attitude  Maneuvers 


C.K.  Carrington*  and  J.L.  Junkinst 
Virginia  Polytechnic  Institute  and  Stale  University,  Blacksburg,  Virginia 


Polynomial  feedback  coatruto  for  hrtHUilt,  aoallnear  spacecraft  aldtwU  ataocoven  art  developed.  A  five- 
body  coafigaratioa  eoasbllog  of  aa  asymmetric  spacecraft  aad  four  reaction  wheels  Is  considered.  Attention  Is 
restricted  to  I  be  momealam  transfer  dass  of  Interaal  control  torques;  this,  la  conjunction  with  the  choice  of 
Eater  parameters  as  attitude  coordinates,  permits  several  Important  order  redact  ion  simplifications.  Three 
namerical  examples  are  Included  to  blast  rate  applications  of  the  concepts  presented. 
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Introduction 

RAPID  large-angle  attitude  maneuvers  have  become  in¬ 
creasingly  important  to  the  success  of  many  current  and 
future  spacecraft  missions.  These  maneuvers  are  characterized 
by  nonlinear  behavior,  however,  resulting  in  a  control  prob¬ 
lem  that  is  likewise  nonlinear.  One  approach  to  feedback  con¬ 
trol  of  nonlinear  motion  is  “gain  scheduling”  in  which  the 
control  history  is  divided  into  segments,  each  determined  by 
its  own  set  of  linear  gains.  A  more  attractive  approach  is  con¬ 
trol  of  the  entire  nonlinear  maneuver  by  a  single  set  of  gains. 

For  the  latter  approach,  a  method  is  presented  whereby  the 
optimal  nonlinear  control  problem  is  solved  in  polynomial 
feedback  form  and  a  suboptimal  control  law  is  determined  by 
truncation.  Currently  there  are  two  approaches  used  to  deter¬ 
mine  the  polynomial  coefficients  for  the  control.  One  is  to  ex¬ 
pand  the  coast-to-go  functional  as  a  polynomial  in  the  states 
and  then  recursively  solve  the  Hamilton-Jacobi-Bellman  equa¬ 
tion,  as  discussed  by  Willenstein,1  Dabbous  and  Ahmed,1  and 
Dwyer  and  Sena.1  In  the  method  used  here,4  the  control  itself 
is  expanded  as  a  polynomial  and  the  coefficients  determined 
recursively  from  the  costate  equations. 

General  Formulation 

Polynomial  state  equations  may  be  written  in  indicia!  nota¬ 
tion  as 

xlmal/Xj  +  cilkxJxk  +  ...  +  biiUj  (i=I,2 . n)  (I) 

where  x,  are  the  states  and  u,  the  controls  to  be  determined. 
Consider  the  optimal  control  problem  of  finding  a  feedback 
control  law  that  brings  the  states  to  zero  while  minimizing  a 
quadratic  performance  index 

1  f '/ 

Jm-T-\  +  riiu,uj  ldl  (2) 

t  o 


Presented  as  Paper  13-2230  at  the  AIAA  Guidance  and  Control 
Conference,  Gatlinburg,  TN,  Aug.  13-17,  1983;  submitted  Sept.  30, 
1983;  revision  submitted  June  18,  1983.  Copyright  ©  American  In¬ 
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The  Hamiltonian  for  this  system  is 

3C  -  Vi  ( q^XtXj  +  r^UjUj  1  +  X,x, 


(3) 


where  it  is  understood  that  x,  is  symbolic  for  the  right-hand 
side  of  Eq.  (1).  The  necessary  conditions  for  a  minimum  pro¬ 
vide  the  state  equation,  Eq.  (1), 


dJC 

ax, 


and  the  costate  equation 


X,-- 


33C 

dx, 


For  unbounded  control 


which  implies 


u,  r y  bkj'Kk 


(4) 


(5) 


(6) 


(7) 


where  r~'  represents  the  elements  of  the  matrix  inverse  of  ru. 
By  assuming  the  costates  can  be  expressed  as  a  polynomial  in 
the  states,  as  in  Ref.  4, 


X,  k  yXj  +  dykXjXt  + .. . 


(8) 


a  nonlinear  feedback  controt  law  is  determined  in  which  kv(t) 
and  *v»  (f)  are  the  control  gains  sought.  By  substituting  Eq. 
(8)  into  Eq.  (S)  and  carrying  out  the  ensuing  algebra,  we  are 
led  to  fl  homogeneous  polynomial  equations  of  the  form 


where 


Iajx;  +  f  J3  J  xr,x*  +...  =0  (9) 

lari  =  function  ( A,B,Q,R.K,K ) 

10)  -  function  ( A,B,C,R,K,D,D ) 


and  K  and  D  are  arrays  whose  elements  are  the  gains  ku  and 
dy,.  Since  Eq.  (9)  must  hold  at  every  point  in  the  state  space,  it 
is  concluded  that  the  functions  in  brackets  must  vanish  in- 


dependency,  so  wc  obtain 


[a(A,B,Q.R,k,K))  =0  (10) 

\&(A,B.C,R.K.b,D)\  =0  (II) 

Equation  (10)  is  a  matrix  differential  equation  determining  the 
linear  feedback  gains;  upon  carrying  through  the  details,  we 
find  that  the  scalar  equations  of  Eq.  (10)  are  precisely  the 
elements  of  the  matrix  Riccati  equation  which  generates  the 
optimal  feedback  control  if  all  nonlinear  terms  in  the  state 
equation  are  absent.  The  solution  for  the  matrix  Riccati  equa¬ 
tion  can  be  determined  by  Potter*  or  Turner’s6  method,  in 
which  an  associated  eigenvalue  problem  is  solved  and  matrix 
exponentials  are  used. 

The  quadratic  feedback  gains  are  determined  by  Eq.  (II), 
which  can  be  rearranged  into  a  set  of  linear  differential  equa¬ 
tions  of  the  form 


=  +  ty] 

where 

[i»]  =  function  ( A,B,C,R,K ) 
(>1  »  function  ( A,B,C,R,K ) 


02) 

03) 


Upon  solving  the  Riccati  eqaution  for  the  linear  gains  (t), 
Eq.  (12)  provides  nonautonomous,  nonhomogeneous,  but 
linear  equations  that  determine  the  quadratic  gains  dljk(t). 
For  the  steady-state  case,  Eqs.  (10)  and  (11)  can  be  solved 
algebraically  for  the  constant  feedback  gains  subject  to 

ka„=dtmk  =  0  (14) 


The  kc„  are  solutions  of  the  algebraic  Riccati  equation,  and 
the  d^  are  obtained  by  setting  d„k  =0  in  Eq.  (12)  and  invert¬ 
ing  the  linear  algebraic  system.  In  the  numerical  examples  con¬ 
sidered  herein,  attention  is  restricted  to  the  constant  gain  case. 


Since  for  n  states  there  are  «J(n+  l)/2  equations  in  Eq. 
(12),  we  do  not  include  the  algebra  of  the  system  considered 
here.  The  emphasis  of  this  discussion  is  the  following 
generalization;  After  solving  the  Riccati  equation  for  the 
linear  gains,  one  is  led  to  sets  of  linear  differential  equations, 
of  the  functional  form  shown  in  Eq.  (12),  that  can  be  solved 
sequentially  to  obtain  the  quadratic  gains,  the  cubic  gains,  and 
so  on,  up  to  any  desired  order.  The  differential  equations  for 
the  gains  of  each  order  depend  upon  the  lower  order  gains. 

Scalar  Example 

Consider  the  optimal  control  problem  of  minimizing  the 
following  performance  index: 

,=-l[V  +  „’,d,  (15) 

2  J<o 


subject  ot  the  state  equation 


Jf=  —  X+M^  +  W 

(16) 

The  costate  equation  is 

X=  -x+X-2<Xjr 

(17) 

and  the  control  is 

u  =  —X 

(18) 

Table  1  Scalar  example 
performance  indices 


Feedback 

order 

Performance 

index 

1 

4314.8 

2 

3796.4 

3 

3725.5 

4 

3717.2 

5 

3716.8 

Fig.  I  Scalar  example;  I,  linear  feedback;  2,  linear  plus  quadratic  *> 

feedback;  },  linear  through  cubic  feedback.  Fig.  i  NASA  standard  four  reaction  wheel  attitude  control  system. 
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Assuming  the  costate  as  a  polynomial  in  state  x. 


\  =  k,x  +  k2x2  +ksx'  -t  ...  (19) 

vi 


i 

•v, 

I 


then  the  coefficient  differential  equations  corresponding  to 
Eqs.  (10)  and  (II)  and  higher  order  terms  are 

Ar,  -  2A,  +  I  -  Af  =0 

*,-3(*1  +  l)*j»-3*,t 

k}-4(kt  +  \)k}  =  2kl-4kit 

*4-5(*,  +  l)*4  =  3*2it,-5*j< 

*'s  -  6(*,  +  I  )ks  -6k1k<  +  3 k]  -  6ktt 


Substitution  of  the  solution  for  z,(r)  into  Eqs.  (21)  and  then 
Eq.  (19)  yields  a  polynomial  feedback  control  law  with  time- 
dependent  coefficients. 

A  numerical  example  for  the  scalar  case  is  included,  using 
t=0.0l  and  tf=  5  s.  The  performance  indices  are  given  in 
Table  1,  and  the  state  variable  and  control  histories  are  given 
in  Fig.  I.  Curves  corresponding  to  fourth-  and  fifth-order 
feedback  are  coincident  with  third  order  and,  hence,  are  not 
plotted.  Fifth-order  polynomial  feedback  essentially  has  con¬ 
verged  to  the  optimal  control  for  this  scalar  problem. 

A  system  for  attitude  control  of  a  spacecraft  with  four  reac¬ 
tion  wheels  is  now  examined. 

Spacecraft  Orientation 

Euler  Parameters 

A  spacecraft  body-fixed  reference  frame  1 6 1  is  related  to  an 
inertial  frame  |ii)  by  the  direction  cosine  matrix  (C(0) ] , 


1 


*«-(«  +  I)  (*i  +  !)*„=(»  +  I)|  kjk„. ,  +  k>kn_  j  + . .. 

+  k*nknn+ 1  ~k„- ,f  I  for  n  even 
k„  -  (n  +1 )  (*,  +  1  )k„  =  (n  +  1 )  ( Ar2Ar„_ ,  +  + ... 


|6|  =  (C(0)](n)  (23) 

where  [  C(0)  ]  is  defined  in  terms  of  the  four  Euler  parameters 
(0O,0,, 02,03).*  These  attitude  variables  are  related  to  the 
body-frame  components  of  the  spacecraft  angular  velocity  u 
by  the  following  kinematic  differential  equations: 


K  +*<«-i>/t*<«*j)/j  + I  fornodd  (20) 

**  Making  the  change  of  variable  from  time  /  to  time-to-go 
r  =  tf—t  and  assuming  a  solution  of  the  form 


k\  —  A|ss  +  Z|  1  , 

*2  =  2fSZl 
*j=*r42j 

(21) 


¥ 


t  where  klss  is  the  steady-state  solution  for  kt,  we  obtain  the 
following  equations: 

I  2(1 +Aiss)z,  -  1  =0 


■~j~  3(1  +  fisstej  =  3(*,ssz}  +  Z|)r 


■4  — 4(1  +Alss)Zj  =4z,z2r  —  2Z|  Jz| 


E 


Y  dz, 

dr 


—  («  +  IKl  +*iss)z»  =  («+  l)IZ|Zn-i«-zfJ  (z2z„.| 


+  Z)Z..2  +  ...+z„/2*„)l  for  n  even 


—  (3|U>|  +  0jU>2  +  03ttfj),  0,  —  '/j  (0O<*>1  —  0)«2  +  02u,j) 

02  =  >//2  (^5<JI  +  0OU>2  ~  0l‘»j),  0j  =  —  Vl  (02U)|  —  (S | CUj  —  0o<i) }) 


with 


(24) 


Momentum  Reference  Frame 

In  addition  to  using  an  arbitrary,  general  inertial  frame 
(til,  a  special  inertial  angular  momentum  frame  |h|  is  in¬ 
troduced  where  is  aligned  with  the  system  angular  momen¬ 
tum  H,  as  discussed  by  Vadali  and  Junkins*  and  Kraige  and 
Junkins.10  The  other  two  unit  vectors  can  be  defined  by  the 
directions  n,  and  ns  assume  after  «i2  is  rotated  to  coincide  with 
H  (see  Fig.  2).  This  reference  frame  can  be  considered  inertial 
if  the  externa!  torques  are  negligible  during  the  maneuver  and 
only  internal  torques  are  present.  As  is  shown  below,  in¬ 
troducing  this  frame  allows  a  use  of  the  angular  momentum 
integral  to  reduce  the  number  of  state  variables. 

The  orientation  of  |6|  with  respect  to  the  momentum 
frame  I A I  is  given  by  the  projection 

(*1  =  [C(6)]|A1  (25) 

where  the  3x3  direction  cosine  matrix  [C)  is  a  function  of 
four  variable  Euler  parameters  (60.6,  ,62 ) .  The  inertial 
frame  |ri)  is  projected  onto  |A|  by 

I  »*  I  =  ( C(or)  ]  ( if  I  (26) 

where  (a0,a,  ,<*2.03 )  are  constant  Euler  parameters  since  both 
|«1  and  |A|  are  inertial.  Using  the  inertial  frame  | m |  com¬ 
ponents  of  the  system  angular  momentum 


dz« 


Cj  -Jj£ - («+  1)0  +  *iss)*»  “  («+  I)|ZIZB-'C-Zf  :(Z2Zr,-t 


+  z,z„-2  + ...  +  z,„-i,/2*(..  n/2  +  t,/2>lfor  n  odd 


(22) 


£s  The  variable  changes  of  Eqs.  (21)  are  generalizations  of  and 
•ere  motivated  by  Refs.  6  and  7. 

Equations  (22)  are  easily  solved,  subject  to  specification  of 
|*he  boundary  conditions;  e.g.,  Ar,  (r)  =  (r)  =  0  at  r  =  0. 


+HK2n1  +  HKini 

the  constant  a,  Euler  parameters  can  be  defined  as 


(27) 


ao=H/H 

L  2H(HL  +//;,)  J 
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The  reaction  wheel  equations  of  motion  are 


Qj  =  0 


«) 


-H  r  H~H*  I  '1 

•' l  2//<//jS, +#/*,)] 


(28) 


The  6,  Euler  parameters  are  then  related  to  the  a,  parameters 
by  the  bilinear,  orthogonal  equation 


'V 

“0 

-“l 

—  a2 

-a," 

00 

►  — 

a. 

“0 

-aj 

«2 

* 

0. 

*2 

“2 

«3 

“0 

~  a1 

02 

.«3, 

.“3 

~«2 

«t 

“0. 

.03, 

(29) 


and  are  related  to  the  body-frame  components  of  u  by  the  dif¬ 
ferential  equations 


- V\  (6j  Wj  +  63023  +  63U3) 

6,  =  'A  (60oj,  -63023  +62“}) 


63  —  ^  (630)3  *4"  6qO)2  63403) 


63  =- ,/i(52ui -6|Uj -6,30)3)  (20) 

It  can  be  shown  from  Eq.  (25)  and  using  the  algebraic  expres¬ 
sions  for  [c(5)  ]  from  Ref.  8  or  10  that 

=  2(6,62  +  606,)*,  +  (*g -6?  +  -5l)62 

+  2(6253  -6063)63  .  (31) 

so  the  body-frame  components  of  the  system  angular  momen¬ 
tum  can  now  be  written  from  H-Hh1  as 


//,  =2(6,63  +  6063  )H 
=  + 


H,  =  2(6363-606,)//  (32) 

Thus  we  have  an  explicit  relationship  to  eliminate  H,  in  terms 
of  6,. 


Spacecraft  and  Reaction  Wheel  Dynamics 

An  arbitrary  asymmetric  spacecraft  with  four  reaction 
wheels  in  NASA  standard  configuration  is  considered  (see  Fig. 
3).  The  system  angular  momentum  H  is  the  sum  of  the 
spacecraft  and  wheel  angular  momenta;  the  body  components 
of  H  and  u  are  related  by 


//=[/•]<-+ (C]r(^]0  (33) 

where  [/*)  is  the  system  inertia  matrix  with  respect  to  the 
body  frame  1 61,  G  a  vector  of  the  four-wheel  angular 
velocities,  and  [J]  the  wheel  axial  moment  of  inertia  matrix 
defined  by  [,/)=diag  \J„ |,  /=  1,2,3, 4.  [C]  is  a  4x3  matrix 
whose  rows  are  the  three  orthogonal  body-frame  components 
of  four  unit  vectors  along  the  wheel  spin  axes. 

Assuming  negligible  external  torques,  the  time  rate  of 
change  of  the  angular  momentum  is  aero,  and  thus  the 
Eulerian  equation  of  motion  is  obtained 

«  =  [/*]<-+  [<?]7’[y)ft+  (<5]tf  =  0  (34) 


where 


M  = 


0  —  <J, 


«j 

-w2 


“2 


0  — u)| 


(35) 


l./)fi+[./]l£]u)=M  (36) 

where  the  four  components  of  the  control  vector  u  are  the  ax¬ 
ial  torques  applied  to  their  respective  wheels  by  the  motors.  To 
eliminate  the  wheel  angular  velocities  G  from  the  spacecraft 
equations  of  motion,  Eq.  (36)  is  multiplied  by  {£]  7  and  then 
substituted  into  Eq.  (34)  to  obtain 

<i«  -IGJU01//+  [<?lr«|  (37) 

where  [CJ  =  (/*  -  Crj£]  ■*,  a  constant  matrix.  Implicit  in 
Eq.  (37)  is  the  substitution  of  Eqs.  (32);  thus,  w=/(w,6,m). 
Notice  the  pleasing  truth  that  rest-to-rest  maneuvers 
(characterized  by  H  =  0)  remove  all  of  the  gyroscopic  terms  in 
Eq.  (37).  It  is  evident  that,  for  this  class  of  maneuvers,  angular 
velocity  control  is  near  trivial  and  attitude  control  is  nonlinear 
only  because  of  kinematic  nonlinearities. 

The  three  equations  of  motion  in  Eq.  (37)  and  the  four  at¬ 
titude  equations  in  Eq.  (30)  will  be  used  to  determine  the  state 
equations. 


Optimal  Feedback  Control  Formulation 
State  Eqnations 

To  obtain  state  equations  of  the  form 

x  =  Ax  +  F(x)+Bn  (38) 

in  which  A  is  a  constant  coefficient  matrix  and  F(x)  a  vector 
function  containing  the  nonlinear  terms,  let  the  state  variables 
be  the  spacecraft  angular  velocities  u,  and  Euler  parameter 
differences  E,\  thus,  the  seven-element  state  vector  is 

x=  (u,  u2  0)3  60  Et  S2  £j  Jr  (39) 

where 

Z,=b,-b,(tj)  (1  =  0, 1,2, 3)  (40) 

These  new  state  variables  introduce  linear  terms  into  the 
dynamic  and  kinematic  equations,  Eqs.  (37)  and  (30), 
respectively. 

The  elements  of  the  A  matrix  for  the  linear  part  of  the  state 
equations  are  found  to  be 

«t.  =812^3  ~8i3«2 

°ij 

<»22=S23«>-*2|W} 

«3I  =832  ^3-833^? 

°33 

0«2  =  -62(r/)/2 

a5i  =2o(,/)/2 

aa  =  (t/)/2 

ati  =^o((/)/2 

=  -M'/V2 

073  =*  60(f/)/2 

«w=0(»=l . 7; 


fli2  ~ 

02l5=S22W3-*23«2 
fl23  =  821^2  "822^? 

0<i  =  —  6,  (f/)/2 
a,  3  =  -  63  (if)/2 
Os2  =  ~  6j  (tj)/2 
°»i  =  «,('/)/2 
**3  —  ~ 6,  (tj)/2 
072  =  ^ I  (tf)/, 2 

7  =  4...  ,7)  (41) 


where  H°  is  Eq.  (32)  evaluated  at  6,(r7)  and  f„  are  the 
elements  of  matrix  | C }  defined  after  Eq.  (37).  Notice  that  au 
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Fig.  4  Case  1:  Spacecraft  angular  velocities. 

are  explicit  functions  of  the  specific  terminal  state  and 
the  magnitude  of  the  system  angular  momentum  H.  B  is  a 
7x4  matrix 

r 


and  the  vector  F(x)  contains  quadratic  and  cubic  terms  in  x,. 

Ptrfonaaace  lades 

Two  quadratic  performance  indices  are  considered 

J.  - -4- \xTQx  +  uTRu  |dt  (43) 

2  J  *q 

and 

J1—lr\,/\xTQx  +  mTWm]At  (44) 

2  Jim 


TIME  (SEC) 

Fig.  5  Case  1:  Wheel  angular  velocities:  t,  linear  feedback;  2,  linear 
plus  quadratic  feedback. 


is  a  3  x  I  vector  containing  the  orthogonal  components  of  the 
vector  sum  of  the  motor  torques.  /,  penalizes  the  four  motor 
torques  and  J}  penalizes  their  projections  on  (he  principal 
axes. 

In  the  numerical  examples,  performance  index  /,  is  used  for 
maneuvers  involving  all  four  wheels.  Performance  index  J2 
may  also  be  used  with  four  wheels,  but  the  wheel  torques  are 
not  unique,  as  shown  below.  All  examples  utilizing  three 
wheels  use  Jt,  in  which  the  4x3  matrix  [£]  in  Eq.  (45)  is 
replaced  by  its  3x3  nonzero  submatrix.  Q  is  a  positive 
semidennite  weighting  matrix,  and  R  and  W  are  identity 
matrices  for  the  examples  considered. 

Feedback  Control 

For  the  performance  index  in  Eq.  (43),  the  optimal  control 
is 

u~-R-'BlX  (46) 

where  B ,  is  the  B  matrix  of  Eq.  (42).  The  analogous  develop¬ 
ment  for  the  performance  index  of  Eq.  (44)  uses  the  state 
equations  in  the  form 

i~Ax  +  F(x)  +  Bj»m  (47) 
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Fif.  •  Case  1:  Control  torques. 


where  B2  is  a  7  x  3  matrix 


■-fer 

0 


The  optimal  control  m  for  J2  is 

m=-W-'Bl\  (49) 

and  the  wheel  torques  «  are  obtained  by  inverting  Eq.  (45). 
This  solution  for  « is  not  unique  unless  [C]  is  a  square  matrix. 


Numerical  Examples 

Several  examples  are  considered  for  an  asymmetric 
spacecraft  with  four  reaction  wheels,  as  shown  in  Fig.  3.  The 
moments  of  inertia  of  the  spacecraft  without  the  wheels  and 
the  wheel  axial  moment  of  inertia  are  given  in  Table  2;  the  in¬ 
ertia  matrix  (/*]  and  the  wheel  geometry  matrix  [£]  are 
given  in  Table  3.  All  examples  were  performed  using  linear 
and  quadratic  steady-state  gains  with  free  Final  time  tf.  The  ex¬ 
amples  start  with  zero  initial  conditions  for  wheel  and 
spacecraft  velocities  and  end  with  large  initial  conditions. 


A  four-wheel  maneuver  using  J ,  and  several  three-wheel 
maneuvers  using  J2  are  executed,  all  with  zero  initial  wheel 
speeds.  Boundary  conditions  are  given  in  Table  4  using  the 
3-1-3  Euler  angles.  Table  5  contains  the  performance  indices. 


Table  2  Moment*  of  lamia, 
bg-m* 

h  *6.215 

h  *5.070 

I,  113.565 

J.  0.05 


Table  3  Inertia  and  wheel  geometry  ms 
*7.212  -  0.2237  -  0.2237 

l/*J  -  -0.2237  *6.067  -  0.2237 

-0.2237  -  0.2237  114.562 

I  0  0 

0  I  0 

|£l« 

0  .0  I 

VI/3  VI/3  VJ/3 


Table  4  Case  1:  Boundary  conditions 


Initial  states 

0.0001 
0.0001 
0.0001 
-m/7. 
-m/3 
-  e/4 
-0.54611 
0.47921 
0.676*7 
0.11820 
-0.33141 
0.46194 
-0.19134 
0.80010 
0.0 
0.0 
0.0 
0.0 


‘Specific  final  boundary  condition*  for  need  not  be 

formally  enforced;  these  are  determined  implicitly  because 
angular  momentum  is  conserved;  i.c.,  for  H  =  const  and 
u(iy)  specified,  0(<y)  is  implicitly  constrained  by  Eq.  (33). 


Table  5  Case  I:  Performance  Indices 


Final  states 

0.0 

0.0 

0.0 

m/2 

m/3 

m/A 

-0.30257 

-0.13976 

0.81747 

046974 

0.33141 

0.46194 

0.19134 

0.80010 


Linear 

feedback 

Linear  plus 
quadratic 
feedback 

Four  wheels, 

Q  i 

6.15126 

6.13691 

ft 

5.76886 

5.62314 

Skew  wheel  off,  J2 

0, 

6.48600 

6.47143 

Oj 

5.92983 

5.73828 

Third  wheel  off,  J2 

Qt 

5.V3377 

5.76176 

Second  wheel  off,  J j 

Qi 

5.92980 

5.76077 

First  wheel  off,  12 

Qi 

5.92962 

5.76071 

Fig.  7  Case  2:  Control  torques:  1,  linear  feedback;  2,  linear  pins 
quadratic  feedback. 


which  are  evaluated  at  120  s  for  the  four-wheel  maneuver 


and  r,-240  s  for  the  three-wheel  maneuvers. 


A  comparison  of  performance  indices  with  full  and  partial 
weighting  on  the  Euler  parameters  is  made.  Q,  indicates  full 
state  penalties  in  the  performance  index,  and  Q2  implies  that 
no  penalties  are  put  on  i0.  Figures  4  and  5  show  the  spacecraft 
and  wheel  angular  velocities  for  the  case  of  three  orthogonal 
wheels  (skew  wheel  off)  and  Q}  in  performance  index  Jlt 
which  produced  the  lowest  three-wheel  performance  index. 
The  control  torques  are  shown  in  Fig.  6. 


Cast  2 


A  rest-to-rest  maneuver  with  nonzero  initial  wheel  speeds  is 
performed  with  three  orthogonal  wheels  (skew  wheel  off)  and 
no  penalty  on  60  in  performance  index  J2.  This  weighting  pro¬ 
duced  the  lowest  performance  index  in  cases  1  and  2.  The 
boundary  conditions  are  given  in  Table  6  and  the  performance 
indices  evaluated  at  tj  -  240 1  are  given  in  Table  7.  Since  the  in¬ 
itial  and  final  Euler  angles  are  the  same  as  in  case  f ,  the  Euler 
parameters  0,  are  the  same,  but  the  system  angular  momen¬ 
tum  is  larger  and,  hence,  the  Euler  parameters  6t  are  different. 
Figures  7-9  plot  the  spacecraft  and  wheel  angular  velocities 
and  control  torques.  Note  that  the  final  wheel  speeds  are  75, 
50,  and  100  rad/s  for  0,,  Q,.  and  Q„  respectively. 


Cast} 

A  three-dimensional  maneuver  is  demonstrated  for  the 
three-wheel  configuration  with  large  initial  spacecraft  angular 
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TIME  (SEC I 

Fig.  I  Case  2:  Spacecraft  aagalar  velocities. 


velocities  and  the  initial  wheel  speeds  of  case  2.  The  boundary 
conditions  for  this  maneuver  are  given  in  Table  8.  Perfor¬ 
mance  index  J2  was  used,  once  with  equal  weights  on  all  Euler 
parameters  5,  and  once  with  no  weight  on  60.  Only  linear  feed¬ 
back  was  used  for  both  performance  indices  to  demonstrate 
the  differences  between  full  and  partial  weighting,  although 
quadratic  feedback  improves  control  performance  when  stable 
linear  control  is  obtained.  Figures  10  and  1 1  give  the  time 
histories  of  the  Euler  parameters,  spacecraft  angular 
velocities,  and  wheel  speeds  for  these  two  cases. 

pbwulcs 

The  performance  indices  in  Tables  5  and  7  were  reduced 
when  quadratic  feedback  was  added  to  the  linear  control.  In 
both  cases  quadratic  feedback  produced  Euler  parameters  that 
reached  their  final  states  earlier  and  resulted  in  a  lower  perfor¬ 
mance  index.  This  also  occurred  with  the  spacecraft  angular 
velocities  in  case  2. 

A  comparison  was  made  between  full  and  partial  Euler 
parameter  weighting  in  the  performance  index,  since  specify¬ 
ing  values  for  three  Euler  parameters  automatically  produces  a 
value  for  the  fourth  when  the  nonlinear  state  equations  are 
used.  No  stability  problems  were  encountered  using  partial 
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WHEEL  angular  velocities 


weighting  as  long  as  spacecraft  velocities  were  small,  but  in 
case  3  partial  weighting  produced  the  tumbling  behavior  in 
Fig.  10.  This  unstable  behavior  occurs  because  the  linear  feed¬ 
back  gains  are  determined  by  only  the  linear  parts  of  the  state 
equations.  Although  the  nonzero  eigenvalues  of 
A  —  BR~'BTK  have  negative  real  parts,  so  that  the  linear 
closed-loop  system  is  stable,  the  nonlinear  terms  in  the  state 
equations  are  large  and  produce  unstable  behavior  in  the 
simulation.  Similar  behavior  is  seen  in  cases  1  and  2,  except 
that  the  quadratic  terms  in  the  state  equations  change  i0  by  a 
smaller  perturbation  that  is  not  destabilizing. 

When  the  departure  motion  of  all  four  Euler  parameters  is 
penalized  in  the  performance  index,  linear  gains  are  deter¬ 
mined  that  will  keep  50  close  to  its  desired  final  state.  The 
quadratic  terms  in  the  state  equations  then  do  not  produce 


destabilizing  corrections  to  the  linear  closed-loop  equations, 
since  bounds  have  been  placed  on  the  states  through  the  per¬ 
formance  index.  As  is  evident  in  comparing  Figs.  10  and  II, 
introducing  the  penalty  on  &0  departure  motion  eliminates  the 
tumbling  motion  and  yields  an  attractive  optimal  maneuver. 

Conclusions 

Polynomial  feedback  on  angular  velocities  and  Euler 
parameters  have  been  used  for  nonlinear  control  of  a 
spacecraft  with  four  reaction  wheels.  A  comparison  of  linear 
and  quadratic  control  was  made,  with  a  reduction  in  the  per¬ 
formance  index  for  quadratic  feedback.  When  using  redun¬ 
dant  attitude  variables,  care  must  be  taken  so  that  the  linear 
gains  (based  upon  linearized  departure  motion  state  equa¬ 
tions)  result  in  modest  violations  of  the  implicit  constraints. 
This  may  be  enforced  by  penalizing  all  states  in  the  perfor¬ 
mance  index  weight  matrix. 
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SPACECRAFT  ATTITUDE  CONTROL  USING 

GENERALIZED  ANGULAR  MOMENTA 

Connie  K.  Carrington  * 

John  L.  Junkins  * 

Suboptimal  nonlinear  state  feedback  control  laws  are  developed  for  space¬ 
craft  attitude  control  using  the  Euler  parameters  and  conjugate  angular 
momenta.  Time-dependent  gains  are  determined  in  closed  form  for  polynomi¬ 
al  feedback  laws,  and  stable  nonlinear  state  feedback  laws  ara  developed 
from  Lyapunov  functions.  The  Lyapunov  laws  are  made  optimal  by  adjusting 
constants  to  minimize  a  quadratic  performance  index.  Numerical  simula¬ 
tions  foi;  single  and  three-axis  large-angle  slew  maneuvers  are  presented. 

INTRODUCTION 

Many  current  and  future  spacecraft  missions  require  rapid  large-angle  reorientation 
maneuvers.  Traditional  feedback  controls  may  not  be  adequate  for  the  accuracy  and 
speed  required  for  on-board,  real-time  implementation,  prompting  investigation  of 
control  laws  using  alternate  state  variables.  Recently,  feedback  control  laws  using 
the  four  Euler  parameters  or  quaternions1,2  and  the  spacecraft  angular  velocities 
have  been  developed  for  attitude  control  (see  Refs.'  3-8). 

Morton9  has  presented  a  new  formulation  of  rigid  body  rotational  dynamics  in  terms 
of  four  generalized  angular  momenta  that  are  conjugate  to  the  Euler  parameters. 
Rotational  motion  is  determined  by  eight  state  equations  that  are  cubic  polynomials 
in  the  states.  The  development  includes  applied  torques,  which  will  be  used  in  this 
paper  for  spacecraft  attitude  control. 

Suboptimal  polynomial  state  feedback  control  laws  8,10  that  minimize  a  quadratic 
performance  index  are  developed  for  rapid  large-angle  spacecraft  maneuvers.  The 
costates  are  written  as  polynomials  in  the  states,  producing  sets  of  differential  equa¬ 
tions  for  the  gains.  These  equations  are  solved  recursively,  since  the  linear  and 
zeroth-order  gains  determine  coefficients  in  the  equations  for  the  quadratic  gains. 
Each  set  of  equations  are  solved  in  closed  form,  producing  gains  that  are  polynomi¬ 
als  in  time.  A  suboptimal  control  law  is  generated  by  truncation  of  the  polynomial 
expansion  in  the  states.  A  numerical  example  using  this  control  law  is  presented 
for  a  single-axis,  large-angle,  spin-down  maneuver. 

Nonlinear  state  feedback  control  laws  utilizing  Lyapunov  functions11  are  also  in¬ 
vestigated  for  single  and  three-axis  maneuvers,  as  in  Refs.  3  and  4.  Unlike  the 
time-dependent  polynomial  control  laws  examined  earlier,  these  feedback  controls 
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University  of  South  Carolina,  Columbia ,  South  Carolina  29208. 

+  Professor  of  Engineering  Science  and  Mechanics ,  Virginia  Polytechnic 
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art*  guaranteed  to  produce  a  closed-loop  system  that  is  asymptotically  stable  in 
the  large.  Constants  in  the  Lyapunov  laws  may  also  be  adjusted  to  minimize  the 
quadratic  performance  index  of  an  optimal  control  formulation.  The  system  re¬ 
sponse  for  the  optimal  Lyapunov  control  law  demonstrates  the  moderate  rise-time, 
short  settling  time,  and  low  overshoot  associated  with  optimally  damped  (f  =  0.707) 
second-order  systems. 


EQUATIONS  OF  MOTION 


The  attitude  control  problem  for  a  rigid  spacecraft  is  governed  by  a  set  of  kinematic 
equations  defining  orientation  of  the  body  with  respect  to  an  inertial  frame,  and 
a  set  of  dynamic  equations  representing  rotational  motion.  The  orientation  of  a 
body-fixed  reference  frame  {b}  to  an  inertial  frame  {n}  is  given  by  the  projection 


m  =  [cm 


where  [C]  is  the  direction  cosine  matrix.  Instead  of  three  Euler  angles,  the  four 
Euler  parameters12  will  be  used  to  parameterize  the  elements  of  [C],  The  Euler 
parameters  are  defined  as 


00  ~  cos(S>/2) 
0,  =  /,  sin($/2) 


where  /,  are  components  of  a  unit  vector  along  the  principal  axis  of  rotation  when 
rotating  from  {n}  to  {6},  and  $  is  the  rotation  angle  for  that  reorientation.  Since 
rotational  motion  has  three  degrees  of  freedom,  the  four  Euler  parameters  are  once 
redundant  and  satisfy  the  constraint 


The  time  derivatives  of  the  Euler  parameters  are  functions  of  the  spacecraft  angular 
velocity  components  w,  and  the  Euler  parameters. 


Rotational  motion  is  governed  by  Euler’s  equations,  which  are  generally  written  in 
terms  of  the  angular  velocity  components  u Consistant  with  the  development  of 
Hamilton’s  canonical  equations12,  Euler’s  equations  are  reformulated  in  terms  of 
generalized  angular  momenta  p,  that  are  conjugate  to  the  Euler  parameters.  They 
are  defined  from  the  rotational  kinetic  energy  T  as  follows 


dT 

Pi  —  •  i  —  0, 1 , 2, 3 

dp, 


The  equations  governing  the  time  derivatives  of  the  Euler  parameters  are  also  re¬ 
formulated  in  terms  of  the  generalized  angular  momenta,  producing  the  following 
eight  equations  of  motion9 
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where  \Q(P) j  and  [7  1  j4  art*  4x4  matrices  defined  as  follows 
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72,  and  J5  are  the  spacecraft  moments  of  inertia  in  a  body-fixed,  principal-axis 
reference  frame.  The  nonzero  elements  of  {u}4  are  the  control  torques  about  the 
principal  axes 

{«}<  =  {  0  «,  «2  «;»  }r  (8) 

Eqs.  (5)  are  Hamilton's  canonical  equations  in  which  {p,  0}  are  the  eight  state 
variables.  These  variables  are  twice  redundant,  with  the  0t  satisfying  Eq.  (3)  and 
the  pf  satisfying  the  following  constraint 


I 


^p,2  =  4//2  (9) 

«-  o 

where  H  is  the  magnitude  of  the  system  angular  momentum  vector.  These  con¬ 
straints  arc  both  integral  properties  of  the  system  equations,  however,  and  hence 
they  do  not  need  to  be  explicitly  enforced  when  defining  the  optimal  control  prob¬ 
lem. 
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OPTIMAL  CONTROL  PROBLEM 
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Formulation 

The  state  equations  are  cubic  polynomials,  so  that  polynomial  feedback  control 
laws  may  be  developed  using  an  optimal  control  formulation8,10.  No  linear  terms 
are  present  in  Eqs.  (5),  however,  and  there  are  fewer  control  variables  than  states,  so 
that  the  algebraic  gain  equations  are  in  many  instances  degenerate.  In  these  cases  no 
solution  can  be  found  for  constant  gains. M  Two  remedies  can  be  considered  for  this 
situation.  The  first  is  to  introduce  linear  terms  into  the  system  equations  by  using 
the  system  errors  as  new  state  variables.  The  new  state  vector  is  z(l)  =  x(/)  -  r(l /), 
where  t(I)  is  the  old  state  vector  ami  r[l j)  is  the  desired  final  state.  Linear  terms 
with  constant  coefficients  r[! j)  appear  in  the  new  state  equations.  The  optimal 
control  problem  is  formulated  using  the  error  vector  c,  and  the  costates  are  written 
as  polynomials  in  c.  Sets  of  equations  for  the*  gains  are  determined  by  equating 
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coefficients  of  like  powers  of  z  in  the  costate  equations.  Constant  linear  gains  can 
be  determined  from  the  algebraic  Riccati  equations  if  at  most  one  zero  pole  occurs 
in  the  linear  part  of  the  state  equations,  and  higher  order  gains  can  be  calculated  re¬ 
cursively  from  the  subsequent  linear  algebraic  equations.  Unfortunately,  the  system 
represented  by  Eqs.  (5)  contain  two  zero  poles,  corresponding  to  0O  and  p0.  The 
eigenvectors  for  these  poles  cannot  be  determined  accurately  enough  to  produce  a 
reasonable  solution  to  the  Riccati  equation  using  Schur’s  method  14,  so  the  constant 
gain  solution  will  not  be  used. 

The  other  solution  technique  for  state  equations  containing  no  linear  terms  is  to  de¬ 
termine  time-dependent  gains  using  the  original  cubic  state  equations.  The  tracking 
problem  can  be  formulated  using  a  performance  index  that  minimizes  the  difference 
between  the  state  and  the  target,  and  a  final  state  penalty  can  be  posed.  The 
costates  are  written  as  polynomials  in  the  original  states  x,  and  the  gain  equations 
are  defined  as  before  from  the  costate  equations.  The  differential  equations  defining 
the  gains  are  simple,  so  that  closed-form  solutions  may  be  found  by  integration. 
The  terminal  boundary  conditions  for  the  costates,  which  are  determined  by  the 
performance  index,  specify  constants  of  integration  in  the  expressions  for  the  gains. 

Given  the  state  equations  in  Eqs.  (5),  the  performance  index  to  be  minimized  is 

J  =  \zT  (t,)Hz(t,)  +  \f'{’?Qz  +  £Ru<)it  (10) 

where  z  is  the  difference  between  the  state  x  and  the  target  state  r 

4(0  =4(0  “£(</)  ,n) 

4(0  =  {  Po  Pi  Pi  PZ  00  01  02  02  } 

The  necessary  conditions  for  optimality  produce  the  following  costate  equations 

;  _  lf^j*(p)K  * 

—  4  dpt  ^J0k  4^ ji\0)lj  ~  9*7 zl 

*  =  -  \  \^r}^  -  V*  -  («) 

i ,  j ,  fc  —  0, 3 
/  =  1 , ...,  8 

where  A,,  and  are  the  costates  corresponding  to  p,  and  0t  respectively;  the  4x4 
matrix  C  is 

cip)  =  [<J(p)][r1l,|Q(p)]r  (13) 

and  the  terms  in  brackets  are  Jacobians.  The  control  is 

a,  =  -2R-'qt(0)\  (14) 

and  the  tern.mal  boundary  conditions  on  the  costates  are 


Assuming  the  costates  are  polynomials  in  the  states. 


{  ^(0  }  =  +  *ti(0*j  +  <*»>*( t)xjxk  +  - 


then  the  boundary  conditions  of  Eq.  (15)  produce  the  following  conditions  for  the 
time-dependent  coefficients  in  the  polynomial  expansion 


$(</)  = 
K(tf)  =  H 


with  all  higher-order  coefficients  going  to  zero  at  t  =  tf. 


By  substituting  Eq.  (16)  into  the  costate  equations  and  equating  coefficients  of  like 
powers  of  z,  we  obtain  the  following  sets  of  equations 


si  —  Qijrj 

Kj  =  -Q,]  +  fiM 

dijk  ~ 


where  f\  and  f2  are  functions  of  lower  order  coefficients.  Eqs.  (18)  can  be  integrated 
in  closed-form,  and  the  constants  of  integration  determined  by  Eqs.  (17).  The  feed¬ 
back  gains  produced  in  this  method  are  polynomials  in  time,  which  are  substituted 
back  into  Eqs.  (16)  and  (14).  A  suboptimal  control  law  is  determined  by  taking  a 
finite  number  of  terms  in  Eq.  (16). 


Single  Axis  Maneuvers 

For  single-axis  maneuvers  about  the  first  principal  axis,  Eqs.  (5)  reduce  to  the 
following  equations 

*i  =  jji*  1*2*4  ”  *2*3)  -  2x4u 

*3  =  Tf  (*i*3*3  -  *1*4)  +  2x3u 

1  ,  (19) 
*3  =  ^(*1*4  ~  *2X3*4) 


*4  =  Jt(X3*3  -*1*3*4) 


where  —  {  Po  p\  0\  }T  and  u  is  the  control  torque  to  be  determined.  For  the 
following  performance  index 


J  -  T{tj)Hz{tf )  +  ^  JJ{ztQz  +  u2}dt 


rJ*~S‘ 


where  z(t)  =  x(t)  -  r{tj)  and  r(lj)  is  the  desired  final  state,  then  the  costate 
equations  are 


Ai  =  -q\j Zj  -  ix2x4  +  A2(x2x3  ~  2x,x4)  +  A3xJ  -  A4X3X4) 


A2  =  —  92j2j  +  4/(^>  (2x2x3  —  X1X4)  -  A2XiX3  -f  A3x3x4  —  A4x3) 


1 


(21) 


A3  =  -93j2j  +  4j(AiX2  -  A2X!X2  +  A3X2X4  +  A4(xjx4  -  2x2x3))  -  2A2u 
1  *  2 


A4  =  -94j2 j  -  ^|(A1x1x2  -  A2x?  +  A3(2xiX4  -  x2x3)  -  A4xtx3)  +  2Atu 
and  the  control  is 


u  =  2(x4A!  -  x3A2). 


(22) 


Aj  and  A2  are  the  costates  corresponding  to  the  conjugate  angular  momenta  po  and 
Pi,  and  A 3  and  A4  correspond  to  Euler  parameters  Pq  and  Pi- 


For  this  example,  the  weighting  matrices  Q  and  H  are 


Q  =  diag{  q ,  q2  93  94  } 
H  =  diag{  h\  h2h$  /i4  } 


The  zeroth-order  equations  are 

si  =  9iri[tf) 

which,  with  the  boundary  conditions  in  Eq.  (17),  have  the  solution 

st(t)  =  - 9.r, (</)(</  -  0  -  htrt(tf) 

The  linear  equations  are 


i  —  9i 

^33  =  “93  +  4a2 


«  =  1,2 


^44  =  “94  +  4sf 


kij  =  0 


These  equations  have  the  solution 


kti  (0  —  9«(*/  —  t)  +  ht  i  —  1,2 

*33(0  =  93 (</  “  0  +  ^3  -  4rl{tj){qlT  +  9aM*/  “  0*  +  ^2 (</  “  0} 

^44(0  =  94(</  “0  +  ^4  “  4r2(tf){q2r  +  9jM*/  “  02  +  h2{t/  -  t)} 


whore 


T  =  (i3f-t:')/3  +  tft2-t2Jt 
23 


(23) 


(24) 


(25) 


(26) 


(27) 


(28) 


Fig.  1  Single- Axis  Conjugate  Angular  Momenta 

1  st  +  ktJXj 

2  4"  ktJXj  -f-  dtji[X ji 
The  quadratic  equations  are 

<^m  =  ^2/2/ 

^123  ~  ~s2 /4/ 
di2i  =  -«i/4/ 

C?134  =  — 4fcji«2  +  *4/4/ 

=  4fcn«i  —  ^3/4/  ^9) 

^213  = 

dj\i  —  ~*i/4/ 

^223  =  al/2/ 

^233  =  4^22*2  "  *4/4/ 

^234  =  ~  4^22*1  +  *3/4/ 

and  dl}it  =  0  where  not  otherwise  specified.  The  terminal  boundary  conditions  for 
these  equations  are  du*(< /)  =  0,  and  the  solutions  contain  linear  through  cubic 
powers  in  time.  Similar  polynomials  in  time  are  found  for  higher  order  gains. 


Fig.  2  Single- Axis  Euler  Parameters 

1  8,  +  kijXj 

2  s,  +  k{jXj  4-  dijkXjXk  . 

Numerical  Example 

A  rapid  large-angle,  spin-down  maneuver  was  simulated  with  I  =  1.00  kg  m2  and 
tf  =  20  sec.  The  boundary  conditions  for  the  maneuver  are  listed  in  Table  1,  and 
the  performance  indices  in  Table  2.  The  performance  index  weights  were 

H  =  Q  =  diag{  0.1  0.1  0.1  0.1  }  (30) 

The  conjugate  angular  momenta  and  Euler  parameters  are  plotted  in  Figs.  1  and 
2.  Fig.  3  contains  the  maneuver  angle  and  control  torque  histories. 

Table  1 

SINGLE- AXIS  BOUNDARY  CONDITIONS 


State 


Initial  State  Final  State 


4>  (rad) 

-ff/3 

ir/3 

u  (rad/sec) 

1.000 

0.000 

0o 

0.866 

0.866 

0i 

-0.500 

0.500 

Po  (kg  m2  rad/sec) 

1.000 

0.000 

Pi  (kg  m2  rad/sec) 

1.732 

0.000 
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Fig.  3  Single- Axis  Maneuver  Angle  and  Control  Torque 

1  4,  +  kijXj 

2  s,  +  kijXj  +  dijkXjXk 

Table  2 

SINGLE-AXIS  PERFORMANCE  INDICES 
POLYNOMIAL  FEEDBACK  CONTROL 


7  H  «  10 


»t  +  kijXj 

Si  +  kijij  -I-  dijkXjX k 


3.13 

0.49 


Note  that  A,  =  Si  +  kijX  j  only  controls  angular  momentum,  and  the  addition  of 
quadratic  terms  are  required  for  attitude  control.  Since  the  control  influence  ma¬ 
trix  contains  the  Euler  parameters,  Aj  =  a,  +  kxjXj + dijkXjXk  results  in  cubic  control 
terms  in  the  state  equations.  This  cubic  feedback  law  provides  a  significant  reduc¬ 
tion  in  the  performance  index  by  reducing  the  Euler  parameters  errors,  even  though 
larger  angular  momentum  values  are  incurred. 


LYAPUNOV  CONTROL  LAWS 
Formulation 

Stable  feedback  control  laws  may  be  determined  for  autonomous  systems  from  Lya¬ 
punov’s  second  method11.  To  apply  this  method,  the  state  equations  are  trans- 
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formed  so  that  the  target,  state  is  the  origin.  This  is  accomplished  by  defining  new 
states  z(t)  =  x(<)  —  r(t f)  as  discussed  in  the  previous  section.  A  scalar  positive  def¬ 
inite  function  V (x)  is  defined,  and  a  control  law  is  found  that  makes  the  total  time 
derivative  of  V (x)  negative  definite.  This  control  law  will  then  be  asymptotically 
stable  and  drive  the  system  to  the  origin. 

Consider  the  positive  definite  function 

V  =  2  £*?  (31) 

r=  1 

2,(0  =  *,-(0  -  >•.(</)  (32) 

X  =  {  P0  Pi  P2  P3  00  01  02  03  }T 

8 

V=A'£z,i,  (33) 

1=1 

xt.  Eqs.  (5)  are  substituted  into  Eq.  (33)  to  give  V  the  following  form 
y  =  /iU)«i  +  h{z)u2  +  /3(z)u3  +  f4{x)  +  /5(x)  +  f6(x)  (34) 

where  uj,  u2,  u3  are  the  control  torques  to  be  determined.  For  asymptotically  stable 
control,  V  must  be  negative  definite.  A  suitable  expression  for  V  is 

v  =  ~{C,ff(z)  +  C2f*(z)  +  CzfUi))  (35) 

where  C,,  C2,  C3  are  constants  to  be  defined  later.  By  equating  Eqs.  (34)  and  (35), 
a  solution  for  the  control  torques  is  found 


The  functions  /j,  /2,  /3  are  quadratic  polynomials  in  s.  and  /«,  /s,  /6  are  cubic 
polynomials  in  z,  so  that  the  quotients  in  Eqs.  (36)  are  well-behaved  as  z  — *  0. 
The  control  laws  are  made  optimal  by  adjusting  the  constants  C i,  C2,  and  C3  to 
minimize  the  quadratic  performance  index  in  Eq.  (10). 

Single  Axis  Maneuvers 

The  large-angle,  spin-down  maneuver  that  was  demonstrated  with  polynomial  feed¬ 
back  control  is  repeated  with  Lyapunov  control.  The  boundary  conditions  are  those 
listed  in  Table  2.  The  Lyapunov  function  is 
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Fig.  4  Single- Axis  Conjugate  Angular  Momenta  and 

Euler  Parameters  Using  Lyapunov  Control 


V  =  8{x3a2  -  xAz\}xi  -  rjyi(x)  -  r2(t/)g2(x)  -  r3{tf)g3[x)  -  rA(tf)gA[x) 


where 


0i(?)  =  {x\X2xA  -  x\x3)II 
g2(x)  =  (x,x2z3  ~  x]xA)/I 
y3(x)  =  (x,xj  -  x2x3xA)/J 
gA[x)  =  [x2x]  -  X\ x3xA) / 1 


To  obtain 


V  =  -C{x3z2  -  xAZi)2  (40) 

let 

This  control  law  produced  the  performance  indices  in  Table  3  for  various  values 
of  C.  These  performance  indices  were  evaluated  at  If  =  30  sec.  The  minimum 
performance  index  occurred  for  C  =  2.9,  and  the  corresponding  system  response  is 
plotted  in  Figs.  4  and  5.  The  conjugate  angular  momenta  and  Euler  parameters  are 
shown  in  Fig.  4.  and  the  maneuver  angle  o  and  control  torque  u  are  plotted  in  Fig. 
5.  The  weighting  matrices  Q  and  //  in  this  example  are  the  identity  matrix. 
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Fig.  5  Single-Axis  Maneuver  Angle  and 
Control  Torque  Using  Lyapunov  Control 

Table  3 

SINGLE- AXIS  PERFORMANCE  INDICES 
LYAPUNOV  CONTROL 

C  Performance  Index 
1.0  7.05 

2.0  3.13 

2.9  2.61 

3.5  2.71 

To  demonstrate  the  variations  in  response  for  changes  in  C,  a  90°  maneuver  was 
executed.  Fig.  6  shows  this  maneuver  angle  for  several  values  of  C.  When  C  is 
smaller  than  its  optimal  value,  the  response  is  similar  to  an  underdamped  second- 
order  system,  and  for  large  values  of  C  the  response  is  overdamped.  For  C  =  1.5,  the 
maneuver  angle  responds  like  an  optimally  damped  (f  =  0.707)  second-order  system, 
and  the  quadratic  performance  index  is  minimized.  Table  4  lists  the  performance 
indices  corresponding  to  the  curves  plotted  in  Fig.  6. 
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Fig.  6  Single- Axis  Maneuver  Angle  Responses 
for  Several  Values  of  C 


Table  4 


SINGLE- AXIS  PERFORMANCE  INDICES 
LYAPUNOV  CONTROL 


Performance  Index 
2.49 
1.59 
1.48 
1.58 
4.70 


By  substituting  Eqs.  (5)  into  the  Lyapunov  function  of  Eq.  (31),  V  becomes 


8{(-Z62i  +  x522  +  X§X3  -  X7*4)uj  + 

(~x7zl  ~  X8Z2  +  *5*3  +  x6zi)u2  + 
(”*i*l  +  X7Z2  ~~  X8Z3  +  *5X4)113  + 

P(‘/)7'l<?(P)][/-,]4|(?(p)]r^  - 
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Fig.  7  Three- Axis  Conjugate  Angular  Momenta  and 
Euler  Parameters  Using  Lyapunov  Control 


A  negative  definite  form  for  V  is 

V  =  -{C^-Xe*]  +  *S*2  +  *823  “  Z7Z4)2  + 

C2{-x7zi  -  x8z2  +  x5z3  +  x6z4)2  +  (43) 

C3(-XgZi  +  X7Z2  -  *623  +  X524)2} 

By  equating  Eqs.  (42)  and  (43),  and  assigning  all  terms  containing  Jj  to  /4(x), 
all  terms  containing  I2  to  fs(x),  and  those  containing  /3  to  /6(x ),  then  the  control 
torques  are  defined  by  Eqs.  (36).  The  functions  f\[z),  f2[z),  and  f3[z)  are  identified 
by  comparing  Eqs.  (35)  and  (43). 
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A  large-angle  three-axis  spjn-down  maneuver  was  executed  using  the  boundary  con¬ 
ditions  listed  in  Table  5.  The  3-1-3  Euler  angles  6,  0  correspond  to  the  given  Euler 
parameters.  Table  6  contains  the  mass  properties  of  the  spacecraft,  and  the  final 
time  is  tf  =  50  sec.  As  in  the  single-axis  maneuvers,  varying  the  constants  C\ ,  C2, 
and  C3  produce  responses  that  exhibit  the  characteristics  of  damped  second-order 
systems.  Several  values  of  these  constants  with  their  corresponding  performance 
indices  are  listed  in  Table  7.  The  minimum  performance  index  occurs  for 

Cl  =  1//; 

c2  =  1/1,  (44) 


C3  =  l//3 

and  Figs.  7  and  8  plot  the  state  variables,  angular  velocities,  and  control  torques 
for  this  case. 
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Fig.  8  Three- Axis  Angular  Velocities  and  Control 
Torques  Using  Lyapunov  Control 

Table  5 

THREE- AXIS  BOUNDARY  CONDITIONS 


State 

Initial  State 

Final  St 

<t>  (rad) 

-tt/2 

rr/2 

0  (rad) 

—  jt/3 

n/3 

tp  (rad) 

—  7t/4 

w/4 

W!  (rad/sec) 

—  .5 

0.0 

u/2  (rad/sec) 

0.3 

0.0 

W3  (rad/sec) 

0.1 

0.0 

00 

-0.331 

0.331 

0i 

0.462 

0.462 

02 

-0.191 

0.191 

02 

0.800 

0.800 

Pc  (kg  m2  rad/sec) 

0.410 

0.000 

Pi  (kg  m2  rad/sec) 

-0.102 

0.000 

P2  (kg  m2  rad/sec) 

-1.050 

0.000 

Pz  (kg  m2  rad/sec) 

-0.022 

0.000 

n  24  2  7  30 


Table  6 


SPACECRAFT  INERTIA 

\  \ 

Axis  Moment  oK  .*ertia  (kg  m2) 

1  1.00 

2  0.83 

3  0.92 

Table  7 

THREE-AXIS  PERFORMANCE  INDICES 


C, 

C2 

c3 

Performance  Index 

1.00 

1.00 

5.00 

4.893 

1.00 

0.83 

0.92 

4.543 

1.00 

1.00 

1.00 

4.368 

1.00 

1.20 

1.09 

4.225 

1.00 

0.60 

0.20 

6.710 

CONCLUSIONS 

State  feedback  control  laws  have  been  formulated  for  a  new  set  of  state  equations 
representing  rigid-body  rotational  motion.  The  polynomial  state  equations  allow 
development  of  polynomial  feedback  controls  using  optimal  control  theory.  The 
absence  of  linear  terms  in  the  state  equations  and  the  double  redundancy  of  the 
state  variables  preclude  the  use  of  constant  gains,  but  time-dependent  gains  are 
determined  in  closed  form.  A  single-axis  example  demonstrates  that  quadratic 
terms  in  the  costates  are  required  for  complete  attitude  control. 

A  second  group  of  control  laws  is  developed  using  Lyapunov  functions.  Although 
nonoptimal,  these  laws  produce  asymptotically  stable  closed-loop  systems,  and  con¬ 
stants  may  be  adjusted  to  minimize  the  quadratic  performance  index  of  an  optimal 
control  formulation.  The  system  response  for  both  single-  and  three-axis  maneuvers 
exhibit  characteristics  of  a  damped  second-order  system;  the  combination  of  con¬ 
stants  that  minimizes  the  performance  index  produces  a  response  associated  with 
an  optimally  damped  (f  =  0.707)  system. 

A  comparision  of  polynomial  feedback  laws  and  Lyapunov  laws  shows  that  Lya¬ 
punov  control  is  easier  to  formulate.  For  systems  involving  more  than  a  modest 
number  of  states,  the  algebra  required  by  polynomial  control  becomes  prohibitive 
and  an  algebraic  manipulator  should  be  used.  There  are  several  advantages  in  us¬ 
ing  the  conjugate  angular  momenta  for  state  variables  instead  of  the  three  angular 
velocities.  Polynomial  feedback  control  using  angular  velocities  was  examined  in 
Refs.  8  and  10,  and  those  state  equations  produced  coefficients  in  the  gain  equa¬ 
tions  that  are  dependent  on  the  final  states.  Although  the  equations  do  not  need  to 
be  reformulated  each  time  the  target  states  were  changed,  a  system  of  simultaneous 
eouations  must  be  resolved  to  find  new  gains.  In  this  paper,  the  time-dependent 


gains  using  the  conjugate  angular  momenta  arc  found  in  closed-form  and  are  quickly 
calculated  f changes  in  the  target  state.  The  cubic  form  of  the  angular  momenta 
equations  also  produce  Lyapunov  control  laws  that  are  well  behaved  as  the  errors 
go  to  zero. 
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ABSTRACT 

A  quasl-analytical  method  is  presented  for  solving  nonlinear,  open-loop, 
optimal  control  problems.  The  approach  combines  a  simple  analytical, 
straightforward  expansion  from  perturbation  methods  with  powerful  numerical 
algorithms  (due  to  Ward  and  Van  Loan)  to  solve  a  series  of  nonhomogeneous, 
linear,  optimal  control  problems.  In  the  past,  the  only  recourse  for  solving 
such  nonlinear  problems  relied  almost  exclusively  on  iterative  numerical 
methods  whereas  the  asymptotic  perturbation  approach  may  produce  accurate 
solutions  to  nonlinear  problems  without  iteration.  The  nonlinear  state  and 
costate  equations  are  derived  from  the  optimal  control  formulation  and 
expanded  in  a  power  series  in  terms  of  a  small  parameter  contained  either 
explicitly  in  the  equations  or  implicitly  in  the  boundary  conditions.  Each 
order  of  the  expansion  is  shown  to  be  governed  by  a  nonhomogeneous,  ordinary 
differential  equation.  Representing  the  generally  non-integrable, 
nonhomogeneous  terms  by  a  finite  Fouries  series,  efficient  matrix  exponential 
algorithms  are  then  used  to  solve  the  system  at  each  order,  where  the  order  of 
the  expansion  is  extended  to  achieve  the  appropriate  precision.  The  asymp¬ 
totic  perturbation  method  is  broadly  applicable  to  weakly  nonlinear  optimal 
control  problems,  including  higher  order  systems  frequently  encountered  in 
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aerospace  vehicle  dynamics  and  control.  A  number  of  numerical  examples 
demonstrating  the  perturbation  approach  are  included. 

OPTIMAL  CONTROL  FORMULATION 

Consider  a  nonlinear  system  of  the  form 

x  +  rx  +  ax  =  Bu  +  f (xtx,u,u,t)  (I) 

where  x  is  an  nxl  vector  of  system  coordinates,  r  is  a  diagonal  matrix  of 
damping  factors,  a  is  a  diagonal  matrix  of  stiffness  factors,  B  is  the  control 
influence  matrix,  u  is  an  mxl  vector  of  controls,  and  f  is  a  vector  containing 
all  nonlinear  terms.  The  form  of  Eq.  (1)  is  based  upon  the  assumption  that 
the  linear  part  of  the  equations  of  motion  have  been  rendered  independent  via 
a  linear  transformation  to  reduce  the  complexity  of  the  optimal  control 
formulation.  If  all  of  the  n  system  coordinates  are  to  be  controlled,  this 
procedure  is  entirely  arbitrary  although  it  will  often  prove  useful  (with 
regard  to  the  validity  of  the  methods  herein).  However,  for  a  control  problem 
in  which  a  subset  of  the  system  linear  modal  coordinates  are  to  be  controlled, 
the  decoupling  procedure  is  necessary.  Furthermore,  it  is  recommended  that 
the  equations  be  nondimensional ized  to  reduce  the  number  of  parameters 
involved  and  to  isolate  the  dimensionless  parameters  critical  to  the  behavior 
of  the  system.  We  shall  proceed  with  the  development  of  the  perturbation 
approach  in  configuration  coordinates,  based  upon  the  assumption  that  the 
linear  part  of  the  equations  of  motion  have  been  decoupled  while  recognizing 
that  the  use  of  dimensionless  coordinates  would  merely  scale  the  coefficient 
matrices. 

Forming  the  state  vector  of  the  system  coordinates,  controls,  and  their 
first  derivatives,  the  state  equation  is 


z  *  Fz  +  OU  + 


where 


z  «  lxT  xT  uT  uT]T 


(2n+2m)  x  L 


I  0  0 
-r  B  0 
0  0  I 
0  0  0 


(2n+2m)  x  (2n+2m) 


(2rv+2m)  x  m 


U  =  u 


m  x  1 


p  *  [0T  fT  0T  0TJT 


(2n+2m)  x  1 


Including  the  controls  and  control  rates  in  the  state  vector  z  allows  us  to 
penalize  large  control  accelerations  in  the  optimal  control  problem.  This 
will  insure  that  the  control  trajectories  generated  will  be  smooth  and 
continuous,  with  prescribed  (usually  zero)  magnitudes  at  the  initiation  and 
completion  of  the  maneuver  (Ref.  1,2).  We  seek  the  optimal  control  trajectory 
that  minimizes  the  quadratic  performance  index 


f 

J  *  A  zTSz I  +  \  [  (zTQz  +  UTRU)dt 
<t=tf  £  t 


where  R  and  S  are  positive  definite,  diagonal  weight  matrices,  and  Q  is  a 
symmetric,  positive  semi-definite  weight  matrix,  where  Q  =  0  is  not 
excluded.  It  is  clear  that  the  matrix  Q  may  be  used  to  weight  not  only  the 
position  and  velocities  of  the  system  coordinates,  but  the  controls  and 
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control  rates  as  well  due  to  the  Inclusion  of  these  variables  in  the  state 
vector  z. 


The  Hamiltonian,  formed  from  the  system  given  by  Eq.  (2)  and  the 
integrand  of  Eq.  (3),  Is  defined  to  be 

H  =  |  (zTQz  +  UTRU)  +  J(Fz  +  DU  +  £)  (4) 

where  the  costates  \  are  a  set  of  undetermined  Lagrange  multipliers. 
Pontryagin's  necessary  conditions  for  determining  the  optimal  control, 
operating  on  the  Hamiltonian,  yield  the  equations 


A  -  —  *  Fz  +  DU  +  £ 


i  -  -  I?  -  -01  -  FTi  -  l^lTi 


0  =  =  RU  +  °Tx  (7) 

with  the  boundary  condition  x_( t^)  =  Sz(t^). 

Solving  Eq.  (7)  for  U  and  substituting  into  Eq.  (5)  reduces  the  optimal 
control  problem  to  two  coupled  first-order,  nonlinear,  ordinary  differential 
equations.  Combining  the  unknowns,  z  and  _x ,  into  a  single  augmented 
state/costate  vector  X,  the  optimal  control  problem  may  be  restated  as 


*  =  AX  +  e { NLT } 


where 


X  .  UT  xT|T 


2(2n+2m)  x  1 


IF  -0R_1DT 
-Q  -FT  . 


2(2n+2m)  x  2(?n+2m) 


{ N  LT }  = 


2(2n+2m)  x  1 


and  where  the  dimensionless  parameter  c  is  a  “bookkeeping"  term  Indicating  the 
numerical  order  of  the  nonlinear  terms. 

Upon  solving  the  Two-Pjint  Boundary-Value  Problem  given  by  Eq.  (8),  the 
state  trajectory,  which  Includes  the  optimal  controls,  may  be  generated  at  any 
point  within  the  time  interval  tQ  <  t  <  tf.  However,  as  a  consequence  of  the 
presence  of  the  nonlinear  terms,  the  system  governed  by  Eq.  (8)  is 
analytically  intractable.  Although  there  are  many  Iterative  methods  available 
for  solving  such  nonlinear  systems,  we  wish  to  construct  an  approach,  using 
the  most  basic  of  the  perturbation  methods,  that  circumvents  the  iterative 
techniques  in  favor  of  a  quasi-analytical  solution  to  the  optimal  control 
problem. 

AN  ASYMPTOTIC  PERTURBATION  METHOD:  THE  PEDESTRIAN  EXPANSION 

For  any  given  weakly  nonlinear,  differential  equation,  it  is  often 
constructive  to  employ  a  straightforward  expansion  from  perturbation  methods 
to  produce  an  approximate  solution  to  the  problem  (Ref.  3).  We  assume  the 
solution  to  Eq.  (8)  may  be  represented  by  a  power  series  in  terms  of  a  small 
dimensionless  parameter  e 

X(t)  =  XQ(t)  +  £ X x ( t )  +  E2X2(t)  +  ...  (9) 

For  small  nonlinearities  (small  e),  the  series  is  expected  to  produce  accurate 
results  where  the  accuracy  will  improve  as  the  nonlinearities,  and 
consequently  the  parameter  c  approach  zero.  Indeed,  in  the  limit,  as  the 
number  of  terms  in  the  series  approaches  infinity,  the  solution  given  by  Eq. 
(9)  will  be  exact  if  the  expansion  is  convergent. 


Substituting  Eq.  (9)  into  Eq.  (8)  the  optimal  control  problem  may  be 


expressed  as 


^  +  eXj  +  t%2  +  0(e3)  *  AX^  +  cAX^  +  e2AX2  +  e{NLT1(XQ)> 
+  £2{NLT2(X0,X1)}  +  0(e3) 


where  the  nonlinear  terms  have  been  expanded  in  a  similar  power  series  and  the 
dependence  of  each  term  upon  the  expansion  variables  (X^)  at  each  order  is 
Indicated.  Equating  terms  with  equivalent  powers  of  e  yields  the  series  of 
equations 


*  *  AX 
— 0  —o 

(ID 

-  AXX  +  {NLT1(X0)> 

(12) 

-2  *  ^2  +  {NLW-1)} 

(13) 

where  for  illustrative  purposes  we  have  included  only  the  equations  through 
2 

order  e  .  However,  we  note  that  the  order  may  be  extended  to  achieve  the 
degree  of  precision  required  for  a  specific  problem.  The  boundary  conditions 
of  the  expansion  variables  are 


W 


(14a) 


M‘o>  ■  j^(t  )}  {”}  1  '>.2.3,...  (14b) 

where  we  recall  that  the  final  conditions  of  the  states  and  costates  are 
related  at  each  order  through  the  boundary  condition  v^(tf)  =  $Ii(tf).  The 
straightforward  expansion  produces  a  strictly  linear  problem  as  the  first 
approximation  (zero  order)  of  the  nonlinear  problem  and  then  provides  a  series 


of  "small  corrections"  (higher-order  terms)  to  account  for  the  effects  of  the 
nonlinearities.  Furthermore,  the  nonhomogeneous  term  in  the  ith  equation  of 
higher  order  (1  =  1,2,...)  is  independent  of  the  expansion  variable  Xi  for 
that  particular  equation.  Upon  solving  Eqs.  (11-13)  sequentially,  the 
significance  of  this  observation  becomes  clear;  the  nonhomogeneous  terms 
constitute  known  functions  of  time.  By  employing  the  perturbation  method,  we 
have,  as  usual,  replaced  an  intractable  nonlinear  problem  with  a  series  of 
nonhomogeneous,  linear,  first-order,  ordinary  differential  equations. 

SOLUTION  OF  THE  NONHOMOGENEOUS  EQUATIONS 

The  solution  of  a  system  of  the  form  given  by  Eqs.  (11-13)  can  be  shown 

to  be 

X  (t)  =  eAt[X  (0)  +  f  e~ATd.(T)dT]  ,  i  =  0,1,2,...  (15) 

1  0 

where 

d^t)  *  0 

d.(x)  =  {NLT . }  ,  i  =  1,2,... 

and  where  tQ  *  0  has  been  assumed  without  loss  of  generality.  Although  Eq. 
(15)  provides  the  solution  for  each  expansion  variable,  evaluating  the 
integral  for  an  arbitrary  function  d^  presents  a  formidable  task  and  may 
require  numerical  integration.  On  the  other  hand,  if  the  nonhomogeneous  term 
may  be  represented  accurately  by  a  continuous  function  of  time,  in  exponential 
form,  the  entire  solution  given  by  Eq.  (15)  may  be  evaluated  using  a  matrix 
exponential . 

Since  we  have  no  closed  form  expressions  for  the  nonhomogeneous  terms,  at 
best  we  can  generate  a  large  set  (k)  of  data  points  which  are  the  "sampled" 
trajectories  of  the  forcing  terms.  A  finite  Fourier  series  of  the  form 
r 

f(t)  =  bn  +  y  a, sin  1w  t  +  b.cos  iu.  t 


(16) 


■fc 


* 


may  then  be  used  to  represent  the  nonhomogeneous  terms  as  a  continuous 
function  of  time.  To  calculate  the  coefficients  of  the  series,  we  use  a  least 
squares  fit  of  the  series  to  the  data  points  describing  the  nonhomogeneous 
function.  It  can  be  shown  that  the  Fourier  series  of  the  jth  element  of  the 
1th  forcing  term  may  be  put  into  the  form 
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(17) 
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where 


“1  0 

1  0 

1 

0 

1  “ 

1  s(t1) 

c(xj)  s (2t 

c(2tj)  ••• 

s(rt1) 

c(rTl) 

A  * 

1  S(t2) 

• 

c(t2)  s(2t2) 

c(2t2) 

• 

s(rt2) 

c(rT2) 

.1  *(Tk) 

c(xk)  s(2tk) 

c(2tk) 

s(r'k> 

c(r,k)_ 

sj  - 

(boj  alj 

blj  a2j  b2j 

...  arj  brj]T 

V 

(d^(°)  d'j(at)  dJ(2at)  . 

..  d^kAt)JT 

T* = 

iui  at  , 

0 

i  =  1,2,3,.. 

•  »k 

at  = 

tf/k 

“o  =  Zr/tf 

c(  ) 

*  cos (  ) 

s(  )  =  sin( 

) 

and  the  notation  d^(kat)  indicates  the  jth  element  of  d^(t)  evaluated  at  t 
=  kat.  The  unknown  coefficients  may  then  be  found  from  the  least  squares 
approximation 

Cj  «  (ATA)-1ATdJ.  (18) 

Note  that  with  appropriate  sample  points  (symmetric  about  t  *  *),  ATA  is 
diagonal  and  the  inverse  is  trivial.  Alternatively,  we  can  develop  these 
coefficients  from  a  discrete  Fourier  transform.  Proceeding  element  by  element 
through  the  ith  forcing  term,  we  can  then  represent  each  element  by  a  finite 
Fourier  series.  It  can  also  be  shown  that  the  forcing  term  may  then  be  given 
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by  the  matrix  exponential  equation 


*i(t>  -  6,e' \ 


where 


Gj  *  l3i  &2  a3  •"  S,'”  •  s  "  2(2lvt2m) 


Sj  *  l»oJ  “0«U  blj  2“o*2J  b2j 
J  *  1,2,3,. ...s 


r“oarJ  brj> 


J 

1! 

“ 


^=(10101...  UT  (2r+l) 


o  =  Block  Diag  (O.Qj.Q-,. ..,Qrl 


(2r+l)  x  (2r+l) 


* 


Hu  )  0 

0 


s,  =  1,2,3,. ..,r 


We  now  have  the  nonhomogeneous  term  represented  by  a  continuous  function 
of  time  given  in  exponential  form.  Consequently,  we  may  use  "Van  Loan's 
identity"  (Ref.  4)  to  produce  the  solution  given  by  Eq.  (15)  using  any 
available  matrix  exponential  algorithms^.  To  accomplish  this,  we  define 


’.■c  :i 

and  Van  Loan  proves  that 
Y  t 

•’  =l  0  «,(t)  . 


where  the  state  transition  sub-matrices  satisfy  the  identities 


♦j(t)  - 

♦„(t)  *  eAt  J1  e'ATG,e°Td 


(22a) 


(22b) 


♦3(t)  -  e 


(22c) 
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Clearly,  Eqs.  (22a,  22b)  may  be  substituted  into  Eq.  (15)  yielding 
*i(t)  =  ♦1(t)Xi(0)  +  ^(t)^ 
where  we  note  that  ^(t)  is  numerically  different  for  each  expansion  variable 
Xj ,  as  indicated  by  the  subscript  (1)  in  Eqs.  (21-23). 


SOLVING  FOR  THE  INITIAL  COSTATES 

Equation  (23)  Is  the  form  of  the  solution  for  each  unknown  variable  (X^), 
however  recalling  the  boundary  conditions  on  the  unknowns,  Eq.  (14),  we 
realize  that  the  Initial  costates,  x^(0),  and  hence  X^(0)  are  as  yet 
undetermined.  The  initial  costates  must  be  found  in  order  for  Eq.  (23)  to 
provide  the  numerical  solution  of  the  optimal  control  problem.  Evaluating  Eq. 
(23)  at  t  =  tf  and  recalling  the  boundary  conditions  from  Eq.  (14),  we  get 

L^f)  )  (1,(0) 


(24) 


It  will  prove  useful  to  write  the  state  transition  matrix  <^(tf)  in 
partitioned  form  and  to  partition  the  last  term  in  Eq.  (24).  Therefore,  we 
define 

■V 


W  = 


Vi^f^  = 


11/*  X 

>1  (tf) 

12 

n 

►ll(*f> 

*22 

*1 

1  ^1 i (tf ^ 

( 

<*2i<V! 

I 

(25a) 


(25b) 


zi(tf)  =  *111(tf)z.(0)  +  *11Z(tf)x.(0)  +ili(tf) 


Substituting  Eq.  (25)  into  Eq.  (24)  yields  the  two  coupled  algebraic  equations 

(26) 

(27) 

in  which  x^(0)  is  tne  only  unknown. 


S^(tf)  =  fZ1(tf )zi (0)  ♦  •Z2(tf)x1(0)  +i2i(tf) 
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Multiplying  Eq.  (26)  by  the  positive  definite  matrix  S,  combining  this 
with  Eq.  (27),  and  collecting  terms  yields 

uf  (tf)  -  s*[2(tf)  11,(0)  =  -  ♦i1(tf)|2,(0) 


+  Sjfe-ii  (tf )  -  ±2i(tf)  (28) 

which  can  easily  be  solved  for  the  Initial  costates  using  any  appropriate 
algorithm.  Now  that  all  of  the  initial  conditions  of  are  known,  Eq.  (23) 
can  be  used  to  produce  the  optimal  control  at  any  time  In  the  interval 
0  <  t  <  tf. 
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RECURSIVE  SOLUTION  OF  THE  STATE  TRAJECTORIES 

Once  the  solution  for  a  given  expansion  variable  is  found,  we  then 
proceed  to  the  next  higher  order.  However,  to  produce  the  Fourier  Series 
approximation  of  the  nonhomogeneous  term  in  the  next  higher  order  requires 
that  we  sample  the  trajectory  of  the  current  order  at  fixed  intervals  of 
time  (At)  throughout  the  maneuver.  Evaluating  the  matrix  exponential 
indicated  in  Eq.  (21)  at  each  time  interval  would  prove  computationally  costly 
and  time  consuming.  An  alternative  procedure  is  to  develop  a  recursive 
formula  for  calculating  the  state  trajectories  whereby  the  matrix  exponential 
is  evaluated  only  once  at  t  =  At.  We  make  use  of  the  exponential  matrix 


recursion 


eA(k+l)at  =  gAatgAkat  (2g) 

to  modify  Eq.  (23).  Applying  the  identity  in  Eq.  (29)  to  the  definitions  in 


Eq.  (21)  produces 

[^((k+ljAt]  [  (k+1) At] 

l  0  $3 ( ( k+1 ) At ! 


'^(At)  *21(at) 

L  o  $3(At) 


(kAt)  $2.(kAt)' 
0  $3(kAt)  . 

(30) 


Carrying  out  the  partitioned  products  indicated  yields  the  three  recursive 
equations  for  the  sub-matrices 


♦  ^  (k+l)at)  -  ♦1(at)«1(kat) 


(31a) 


♦2.[(k*l)Atl  =  *1(at)«21(kat)  +  $21(At)»3(kAt) 


(31b) 


♦3((k+l)at)  *  ♦3(at)«3(kat) 


(31c) 


where  *^(0)  «  I,  *2i(0)  *  °*  and  *3(0)  = 

Similarly,  evaluating  Eq.  (23)  at  t  =  (k+l)At  gives 

^[(k+ljat)  =  «1((k+l)At]X.(0)  +  «2^  l  (k+ljatjc^  (32) 

We  can  then  simplify  Eq.  (32)  by  defining  the  vectors  Vjit  j  =  1,2,3  to  be 

Xjj  l  (k+l)At]  =  «1l(k+l)AtlX.(0) 

v2i((k+l)Atj  =  ®2il(k+l)AtJao 

v3i((k+l)At]  =  4» 3 ( (k+1) Atjc^  (33) 

Substituting  Eq.  (31)  into  Eq.  (33)  and  substituting  this  result  into  Eq.  (32) 
yields  the  recursive  formula 

X.[(k+l)Atl  «  v^Uk+ljAtl  +  v2i((k+l)At)  (34) 


v^Kk+ljAt)  =  *li(At)vli(kAt), 


vn(0)  =  ^(0) 


— 2i [ (k+l)At  1  =  ^(Atjv^fkat)  +  «2i(4t)v3i(kAt),  v2i(0)  =  0 


v3il(k+l)Atl  =  •3(At)v31(kAt), 


^3i (°)  =  90 
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As  a  result  of  Eq.  (34),  the  matrix  exponential  must  be  calculated  only 
twice  for  each  X^;  once  at  t  =  t^  for  solving  the  initial  costates,  and  once 
at  t  *  at  for  the  recursive  formula.  After  completing  the  solutions  of  each 
X.j  in  the  series,  the  optimal  control  trajectory  is  produced  by  combining  the 
trajectories  from  each  expansion  variable  In  the  series  given  by  Eq.  (9).  We 
again  stress  the  fact  that  the  solution  to  the  nonlinear  optimal  control 
problem  has  been  produced  by  solving  a  series  of  strictly  linear,  constant 
coefficient  subproblems  without  the  need  for  iterative  techniques.  We  shall 
illustrate  the  effectiveness  of  the  perturbation  method  with  numerical 
examples  of  low  order  systems. 

NUMERICAL  EXAMPLES 
Case  1: 

To  demonstrate  the  perturbation  method,  the  first  example  of  a  nonlinear 
system  is  a  scalar  problem  with  both  quadratic  and  cubic  nonlinear  terms.  The 
system,  in  configuration  space,  is  given  by 

x  +  cx  +  kx  =  u  +  c(bux  -  ax^)  (35) 

where  c  =  0.1  6  =  1.0  E  =  0.1 

k  =  1.0  a  =  0 . 5 

x(0)  =  1.0  x(Q)  =0  tf  =  2 

The  objective  is  to  determine  the  optimal  controls  that  will  drive  the 
variable  (x)  to  zero  (with  a  final  velocity  of  zero)  in  a  two  second  time 
interval.  To  verify  that  the  system  is  weakly  nonlinear,  the  nondimensional 
form  of  Eq.  (35)  can  be  shown  to  be 

••  _  3 

n  +  6^n  +  ti  =  u  +  S^Un  -  (36) 

with  6^  =  0.1  =  0.1  =  0.05 

where  n  is  the  dimensionless  position  coordinate  and  u  is  the  dimensionless 
control  force.  Clearly,  from  Eq.  (36),  the  system  is  lightly  damped  with  weak 


nonlinearities.  We  shall  choose  to  penalize  only  the  control  accelerations 
and  the  final  state  in  the  performance  measure,  Eq.  (3),  and  as  such  we  let  R 
*  1,  Q  *  0,  and  S  *  10^(I|.  The  effect  of  the  large  weight  matrix  (S)  Is  to 
rigidly  enforce  the  final  conditions  in  the  optimal  control  problem.  In  the 
vernacular  of  optimal  control  theory  this  example  is  a  fixed  time,  fixed  final 
state  optimal  control  problem  (Ref.  6).  Numerically,  we  shall  solve  this 
problem  in  configuration  space,  and  as  such  the  matrix,  F,  is  given  by 


Similarly,  the  vector  of  nonlinear  terms  in  Eq.  (8)  can  be  shown  to  be 


where  x ^  is  the  second  element  of  the  costate  vector  We  evaluate  the 
effectiveness  of  the  optimal  control  approximations  by  integrating  Eq.  (35), 
numerically,  using  a  4-cycle  Runge-Kutta  routine  and  examining  the  final 
boundary  condition  errors  of  the  numerically  integrated  solution. 

A  second  order  expansion  in  the  power  series  yields  a  final  condition 
x(tf)  =  -0.000322  from  the  integrated  equation  of  motion.  While  not  exactly 
zero,  the  error  Is  less  than  0.04%.  By  comparison,  the  linearized  optimal 
control,  obtained  by  dropping  the  nonlinear  terms  (note  that  this  is  also  the 
zeroth  order  expansion  variable)  produces  x(tf)  =  -0.0402  or  an  error  of  over 
4%.  The  perturbation  approach  reduced  the  error  by  two  orders  of  magnitude 


8 


for  a  second-order  exansion,  demonstrating  the  effectiveness  of  the  method. 

The  trajectories  of  the  position  and  control  are  shown  in  Fig.  1,  where  each 
profile  exhibits  the  smooth,  continuous  behavior  expected  of  an  optimal 
solution  of  this  problem. 

Case  2: 

For  the  second  example,  the  perturbation  method  is  applied  to  the  second- 
order  system 


where 


*  klxl  '  U1  *  ‘Vlx2 

+  k^X^  =  +  ec^X^Xp 

1  =  0.1  kj  =  1.0 

al  =  1.0 

2  =  0.1  k2  15  0.5 

=  0.5 

^(0)  =  1.0  i1(0)  =  0 

e  =  0.1 

x2(0)  =  2.0  x2(0)  =  0 

In  nondimensional ,  matrix  format,  the  equations  of  motion  are 

■0.1  01  fl  0  ((O.lKn-  ) 

jn+  n.  +  _n  =  U:  +  <  >  (38) 

.  0  0.1J  l0  0.5 J  f(0.05)n1ri2) 

where  again  we  see  that  Eq.  (38)  is  a  lightly  damped,  weakly  nonlinear 

system.  We  shall  choose  to  penalize  the  final  states  and  the  second 

derivatives  of  the  controls  in  the  performance  index.  Mathematically,  we 

state  this  by  setting  R  =  I,  Q  =  0,  and  S  =  1020 [  I ) .  Proceeding  as  with  Case 

1,  the  results  of  a  second  order  expansior  are  compared  to  the  linearized 

optimal  control  problem  as  shown. 

X)(tf)  x?(tf) 

Linearized  optimal  control  approximation:  0.135  0.0868 

Second  order  control  approximation:  0.0000166  0.00000945 

The  final  position  errors  for  the  zeroth  order  control  approximation  are  13.5% 
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and  4.3%  respectively  for  the  state  variables  xj  and  X2*  The  control 
determined  from  the  second  order  perturbation  expansion  reduces  the  errors  by 
approximately  four  orders  of  magnitude.  Such  explosive  convergence  is  not 
typical  but  it  does  demonstrate  how  well  the  perturbation  method  may  solve 
open  loop  optimal  control  problems. 
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Case  3: 


As  a  final  example,  we  wish  to  test  the  method  with  a  system  containing 
larger  nonlinearities.  To  accomplish  this,  we  shall  use  the  same  system  as 
used  in  Case  2  with  the  following  parameter  changes: 


°1  =  2  xl<°)  =  3  e  =  °*4  a2 

3  x2(0)  =  2 

In  dimensionless  form. 

the  system  is  given  by 

ro.i  01 

•  r1  0 1  .  i 

1(2.4) n, no) 

-  + 

a  +  1  a  =  y.  +  1 

1  1  (39) 

L  0  O.lJ 

Lo  0.5)  | 

[(3. 6)^2  j 

We  notice  immediately  that  this  system  is  strongly  nonlinear,  and  would  not 

expect  controls  from  the  zeroth  order  solution  to  accurately  approximate  the 

actual  optimal  control 

.  However,  we  shall  proceed  with  the  perturbation 

approach  while  recognizing  that  this  is  a  significant  test  of  the  method.  The 

final  conditions  are  shown  in  Table  I  for  controls  computed  from  expansions  of 

zero  (linearized  system)  through  sixth  order. 

TABLE  1  FINAL  STATE  ERRORS 

Approximation  x.(t.)  > 

:2(tf) 

Order 

X  1 

0 

110.00 

168.86 

1 

4.164 

6.859 

2 

0.3654 

0.6124 

3 

-0.01957 

-0.03416 

4 

-0.01001 

-0.01677 

5 

0.0006371 

0.001106 

6 

0.0002444 

0.0003835 

The  controls  from  the  sixth-order  expansion  produces  very  accurate 
results  with  errors  substantially  less  than  0.04%  for  both  coordinates.  The 
position  and  control  profiles  are  shown  in  Figs.  2  and  3  respectively;  each 
trajectory  is  a  smooth  and  continuous  path  to  the  origin.  The  excellent 
convergence  is  achieved  for  this  problem  In  which  the  nonlinear  terms  are  of  a 
significant  magnitude.  In  Ref.  7,  we  document  analogous  results  for  maneuvers 
of  a  flexible  spacecraft;  we  have  achieved  reliable  convergence  for  a  system 
having  order  42. 

CONCLUSIONS 

A  procedure  for  solving  nonlinear,  open  loop,  optimal  control  problems 
has  been  presented.  In  this  approach,  an  asymptotic  perturbation  method  is 
applied,  thereby  obtaining  a  solution  process  without  the  traditional 
dependence  on  iterative  numerical  methods.  The  nonlinear  system  is 
"separated"  into  a  set  of  nonhomogeneous,  linear,  optimal  control  problems 
that  may  be  solved  sequentially.  Upon  combining  the  solutions  of  the 
subproblems  in  a  straightforward  power  series,  an  optimal  control  for  the 
nonlinear  system  is  generated.  This  novel  process  for  solving  nonlinear 
optimal  control  problems  is  a  result  of  the  marriage  of  a  simple  analytical 
technique  (the  perturbation  method)  and  a  powerful  numerical  algorithm  (the 
matrix  exponential). 

Although  the  asymptotic  perturbation  method  was  conceived  as  a  solution 
process  for  weakly  nonlinear  problems,  the  method  has  demonstrated 
extraordinary  effectiveness  when  applied  to  many  strongly  nonlinear  problems 
such  as  the  system  presented  in  Case  3.  Certainly,  the  perturbation  method 
will  not  produce  accurate  results  for  all  nonlinear  systems.  However,  the 
family  of  nonlinear  problems  to  which  the  method  is  effective  is  considerably 
larger  than  we  initially  expected.  We  therefore  anticipate  that  this 


asymptotic  perturbation  method  will  be  found  to  be  broadly  applicable  to  a 
large  family  of  generally  nonlinear  problems,  including  higher  dimensioned 
systems. 
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A  QUASI -ANALYTICAL  METHOD  FOR  COMPUTING 
NONLINEAR  ATTITUDE  MANEUVER  CONTROLS 
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R.C.  Thompson  ,  J.L.  Junkins  ,  and  J.D.  Turner 

A  quasi-analy tical  method  is  presented  for  solving 
nonlinear,  open  loop,  optimal  control  problems. 

Upon  applying  the  most  basic  of  the  asymptotic 
expansions  from  perturbation  methods  (the  pedestrian 
or  power  series  expansion),  the  nonlinear  control 
problem  is  replaced  by  a  sequence  of  linear,  non- 
homogeneous  problems.  In  contrast  to  the  usual 
emphasis  (in  perturbation  methods)  upon  analytical 
solutions  of  the  sequence  of  linear  systems,  we 
show  that  this  sequence  of  problems  can  be  solved 
sequentially  using  efficient  numerical  methods. 

The  nonhomogeneous  terms  are  represented  by  a  finite 
Fourier  series,  allowing  the  use  of  matrix  expo¬ 
nential  algorithms  due  to  Ward  and  Van  Loan  to  be 
used  to  solve  the  system  at  each  order.  In  prin¬ 
ciple,  the  order  of  the  expansion  may  be  extended 
indefinitely,  however  numerical  difficulties  will 
arise  at  some  point.  Historically,  solutions  to 
nonlinear  control  problems  have  relied  almost 
exclusively  on  iteration  methods  (using  one  of 
several  available  algorithms);  however,  the  pertur¬ 
bation  method  presented  here  often  produces  very 
accurate  solutions  to  nonlinear  problems  without 
iteration.  The  perturbation  method  is  broadly 
applicable  to  a  variety  of  optimal  control  problems 
including  large  degree  of  freedom  systems.  Numerical 
examples  of  multi-axis  attitude  maneuvers  of  a  rigid 
spacecraft  are  presented. 
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OPTIMAL  CONTROL  FORMULATION 


Consider  a  system  governed  by  a  nonlinear,  vector 
equation  with  initial  and  terminal  boundary  conditions  of 
the  form 


z_  -  F_z  +  DU_  +  E£ 
z(tQ)  =  a  £_(  t  f )  -  0 


(1) 


where  z_  is  an  nxl  state  vector,  U  is  an  mxl  vector  of 
controls,  e£  is  an  nxl  vector  containing  all  nonlinear 
terms  and  e  is  a  small  dimensionless  parameter.  We  note 
that  many  systems  require  some  algebraic  manipulations 
in  order  to  represent  the  system  in  the  form  of  Eq.  (1). 

For  example,  the  system  may  require  a  transformation  to 
modal  coordinates  in  order  to  produce  an  independent  set 
of  equations.  Furthermore,  the  state  vector  may  be  aug¬ 
mented  to  contain  "pseudo-state"  variables  such  as  the 
control  magnitude  and  control  rate  in  order  to  determine 
more  realistic  or  desirable  control  profiles  and  (in  the 
case  of  flexible  bodies)  reduce  spillover  into  the  higher 
frequency  modes.  (Ref.  1)  If  such  an  approach  is  taken, 
then  the  vector  U  may  not  be  the  applied  control  but  instead 
will  contain  higher  derivatives  of  the  control  variables. 
Because  the  perturbation  method  that  we  present  depends 
upon  a  power  series  expansion,  it  is  important  to  expand 
any  transcendental  functions  into  a  truncated  Taylor's 
series.  Therefore,  if  a  given  system  contains  transcen¬ 
dental  functions,  we  must  start  with  an  approximation  of 
that  system  as  a  consequence  of  this  procedure,  where  of 
course  the  accuracy  of  the  approximation  will  depend  upon 
the  number  of  terms  included  in  the  formulation. 

We  seek  the  control  trajectory  that  minimizes  the 
f unc  t iona 1 

J  =  |£TSz|t  =  t  +  I/VQ£  +  UTRU)dt  (2) 

^  C  o 

subject  to  the  condition  that  Eq.  (1)  is  satisfied;  where 
R  and  S  are  positive  definite,  diagonal  weight  matrices, 
and  Q  is  a  symmetric,  positive  seni-def inite  weight  matrix. 
The  Hamiltonian,  formed  from  the  integrand  of  Eq.  (2)  and 
the  system  equations  is  given  by 


where  the  cos  tales,  designated  as  are  a  set  of  n 
undetermined  Lagrange  multipliers.  Pont  ryagin  '  s  necessary 
conditions  for  determining  the  optimal  control  yield  the 
threeequations 

3H 

Z  ■=  yy  ■  Fz_  +  DIJ  +  E£  (4) 

•  3H  T  3  P  T 

A  -  -  -  -Q*  -  fTA  -e[T|]Ti  (5) 

9  H  T 

0  -  yfj  =  RU  +  nl\  (6) 

and  the  boundary  condition  ^(t^.)  =  S£(t^). 

Solving  Eq  .  (6)  for  JJ  and  substituting  into  Eq  .  (4) 

reduces  the  problem  of  determining  the  optimal  control  to 
solving  two  coupled,  nonlinear,  first-order  ordinary  differ¬ 
ential  equations.  The  two  equations  may  then  be  combined  by 
defining  the  augmented  s t a t e / c o s t a t e  vector,  X,  such  that 


X  =  AX  +  e{XLT  } 


where 


r  T  %T,T 

[z  A  ] 


2n  x  1 


_  1  T 

F  -DR  D 


2n  x  2  n 


{ NLT }  = 


r 

->sTl i' 


2n  x  1 


We  have  now  reduced  the  optimal  control  problem  to  a  Two- 
Point  Boundary  Value  Problem  with  split  boundary  conditions 
on  z_  and  the  condition  that  X/tf)  =  S£(tf).  As  a  consequence 
of  the  existence  of  the  vector  of  nonlinear  terms,  designated 
{NLT},  the  problem  given  by  Eq.  (7)  is  analytically  intract¬ 
able,  and  we  would  ordinarily  rely  upon  iterative  methods  to 


complete  the  solution.  However,  we  wish  to  construct  a 
quasi-analy t  ica 1  approach  that  eliminates  the  reliance  upon 
iteration. 


THE  PEDESTRIAN  EXPANSION 

For  a  weakly  nonlinear  system,  as  indicated  by  the 
"bookkeeping"  term  e,  it  is  constructive  to  apply  a  straight- 
fo  ard  power  series  expansion  to  produce  an  approximate 
solution  to  the  problem.  (Ref.  3)  Let  the  solution  to 
Eq .  ( 7 )  be  given  by 


X(t)  -  X  (t)  +  tX.(t)  +  e 2X _ ( t )  +  ...  (8) 

—  — o  —1  —l 

For  small  non  1 inea r  i  t  i  e s  ,  and  consequently  small  e,  the 
series  will  produce  accurate  results  and  the  accuracy  will 
improve  as  the  non  1 ine ar i t i e s  approach  zero.  Substituting 
Eq.  (8)  into  Eq.  (7)  gives 

XQ  +  eXj  + 1 2X 2  +  0 ( e 3 )  =  AXq  +  eAXj  +  e2AX2  + 

+  e2{ NLT2(Xq,X1) }  +  0 ( e  3) 

where  the  nonlinear  terms,  { NLT  .  }  ,  have  been  expanded  in  a 
similar  power  series  and  the  functional  dependence  of  each 
term  on  the  expansion  variables  (X.)  is  indicated  in  each 
term.  Equating  like  powers  of  e  yields  the  series  of  equa¬ 
tions 


-o  "  A-o 

(10) 

X  =  AXj 

+  { NLT  j  (XQ)  } 

(11) 

X 2  =  AX 2 

+  { nlt2(x0,x1) } 

(12) 

For  illustrative  purposes,  we  have  included  only  the  terms 
through  second  order;  however  the  series  could  be  continued 
to  higher  orders  if  necessary  to  achieve  the  precision  appro¬ 
priate  for  a  specific  problem.  In  Ref  [8],  we  show  several 
examples  wherein  5th  and  6th  order  expansions  were  routinely 
computed;  analytical  perturbation  solutions  above  third  order 
are  very  rare.  The  same  procedure,  when  applied  to  the 
boundary  conditions  yields 


e{  NLT  1  (X0)  } 
(9) 


(  1  3) 


-0 ( t0) 


and  we  recall  that  the  final  conditions  of  the  states  and 
costates  are  related  at  each  order  by  _X  .  ( t  ^  )  -  Sz^^(t^). 

Note  that  the  nonhomogeneous  terms  in  the  ith  equation  of 
Eqs.  (10-12)  are  dependent  only  upon  the  preceeding  i-1 
expansion  variables.  Therefore,  the  nonhomogeneous  term 
in  each  equation  constitutes  a  known  function  of  time  and 
Eqs.  (10-12)  can  be  solved  sequentially.  The  perturbation 
method  has,  as  usual,  replaced  the  original  nonlinear  prob¬ 
lem  with  a  series  of  linear,  nonhomogeneous  problems. 

The  solution  of  Eqs.  (10-12)  can  be  shown  to  be 

X  (t)  =  eAt[X  (0)  +  f  e'ATd.(T)dt]  i-0,1,2,...  (15) 

1  J0  1 

where 

dQ(-r)  =  2 

d  .  (t)  =  { NLT  .  }  1=1,2,...  (16) 

—  l  l 

and  where  we  have  assumed  tn=0  without  loss  of  generality. 
Evaluating  the  integral  in  Eq.  (15)  may  be  accomplished  in 
any  number  of  ways.  Direct  numerical  integration  could  be 
used  but  a  more  attractive  method  would  be  to  use  a  finite 
Fourier  series  approximation  of  the  integrand  (or  some  other 
orthogonal  series)  and  then  integrate  the  series  term  by 
term.  The  latter  method  will  provide  us  with  a  model  of 
the  integral  as  a  continuous  function  of  time,  an  advantage 
when  calculating  the  trajectories  of  the  expansion  variables 
Another  alternative,  and  perhaps  the  most  elegant,  is  to  use 
a  matrix  exponential  to  calculate  the  entire  integral  shown 
in  Eq  .  (15). 

The  matrix  exponential  method  requires  the  nonlinear 
terms,  d.(x),  to  be  represented  by  a  continuous  function  of 
time  in  exponential  form.  Since  we  have  no  closed  form 
expression  for  the  nonlinear  terms,  we  can  at  best  generate 
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a  large  set  (k)  of  data  points  which  are  sampled  values  of 
the  nonhomogeneous  terms.  A  finite  l'our'er  series,  which 
can  always  be  put  in  exponential  form,  is  used  as  the 
continuous  time  representation,  and  is  given  by 

r 

d(t)  =  +  L  a i  sin  iw^t  +  bi  cos  iw^t  (16) 

To  calculate  the  series  coefficients,  we  use  a  least  squares 
fit  of  the  series  to  the  sampled  data  sets.  It  can  be  shown 
that  the  Fourier  series  of  the  jth  element  of  the  ith  non¬ 
linear  term  may  be  given  by 

8c .  =  d  .  (17) 

~3  ~3 

where  the  matrix  8  contains  the  trigonometric  functions 
evaluated  at  the  sample  times,  _c  .  contains  the  unknown 
series  coefficients,  and  ci  .  contains  the  sampled  data  of 
the  jth  element  of  d.  (t),  ils  given  by 

d.  =  [dj(0)  dj(At)  d|  ( 2  A  t )  ...  d|(kAt)]T  (18) 

where  At  =  t,/k.  The  unknown  coefficients  are  then  deter¬ 
mined  from  tne  least  squares  approximation. 

Cj  =  (8T8)"18Tdj  (19) 

Note  that  w^th  appropriate  sample  points  (symmetric  about 
u)QAt  =  tt)  8  8  is  diagonal  and  the  inverse  is  trivial. 
Alternatively,  we  can  develop  these  coefficients  from  a 
discrete  Fourier  transform.  Proceeding  element  by  element 

(j  — 1,2,3 . 2n)  through  the  vector  of  nonlinear  terms, 

each  element  is  represented  by  a  finite  Fourier  series. 

The  collection  of  series  representations  may  then  be  put 
into  the  form 

d.(t)  =  G.ent£()  (20) 

where  the  jth  row  of  G .  contains  the  series  coefficients 
for  the  jth  element  of  the  nonlinear  term,  and  the  matrix, 

Q,  is  defined  to  be 

(1  =  Block  Diag  [0  J 


wli  e  r  e 


i  -  1.2,3,  ..  ,r 


and  Is  a  constant  vector  used  to  "separate"  the  trigono¬ 
metric  functions  calculated  in  Eq .  (20). 

At  this  point,  we  have  a  continuous  time  represen¬ 
tation  of  the  nonlinear  term  in  the  form  of  a  matrix 
exponential.  We  now  employ  Van  Loan's  identity  (Ref.  A) 
which  will  produce  the  integral  in  Eq .  (15)  from  a  matrix 

exponential.  We  shall  define 


Y  .  = 

l 


A  G 


0  0 


it  =  rv*>  *2i(t)' 

.  °  *3(t)  _ 


(21a) 


(21b) 


where  we  note  that  4  (t)  is  different  for  each  X..  The 

definitions  of  the  state  transition  matrices  (4.)1are 
(from  Van  Loan)  1 


$j(t)  =  e 


(22a) 


4>  (t)  =  eAt  f  e  AlG.e^Tdx 

21  J0 


4  3  (  t  )  =  e' 


(22b) 


(22c) 


We  may  then  rewrite  Eq.  (15)  using  the  definitions  from 
Eq.  (22)  to  get 


X.  (t)  =  (t)X.(O)  +  4  2  3 ( t ) £q 


In  order  to  use  Eq.  (23)  to  produce  the  trajectory  of  X., 

we  must  first  determine  the  n  unknown  initial  conditions 

of  X.  (i.e.  the  initial  costates  X.). 

—  l  -  i 


ftt 


SOLVING  THE  INITIAL  COSTATES 

We  have  noted  that  the  costate  boundary  condition 
resulting  from  Pontryagin's  principle  is  _X.(t^)  ■  Sz^(tj-). 
Until  now,  this  boundary  condition  has  not1been  used, 
however,  it  will  be  instrumental  in  completing  the  imme¬ 
diate  task  of  solving  for  the  initial  costates.  Evaluating 
Eq .  (23)  at  t«*t  ^  and  substituting  the  boundary  condition 

into  the  result  gives 


Si^f) 


£  (  tf  ) 

S  z^i  ( tf ) 


z  (0) 

$1  (tf }  '  x  (Q)  +  *2i(tf^0  (24) 


To  simplify  the  notation,  let  us  partition  the  components 
of  Eq.  (24)  as 


*l(tf>  = 


(25a) 


*2i  (  C  f  ^0 


(25b) 


Substituting  Eq .  (25)  into  Eq .  (24)  leads  to  the  two  coupled 

algebraic  equations 

zi(tf)  =  1  (tf)z  .  (0)  +  $,12(tf)X.(°)  +  iu(tf)  (26) 

Sz.(tf)  =  $21(tf)z.(0)  +  $>22(tf  n.(0)  +  i2i(tf)  (27) 

where  X_.(0)  is  the  only  unknown  term.  Finally,  multiplying 
Eq .  (2ot  by  S,  substituting  into  Eq .  (27),  and  collecting 

the  terms  produces 

[»|2(cf)  -  S«]2(tr)  1X^0)  -  [  Sf  ]  '<<:(-)  -  ♦f(tf)U1(0) 


+  'Sli'1!1  '  i2i(tf> 


.■  <■ 


By  using  any  appropriate  algorithm,  we  can  solve  the  linear, 
algebraic  system  given  by  Eq  .  (28).  Substituting  this  solu¬ 
tion  into  Eq  .  (23),  we  now  have  the  solution  for  X.^(t), 

RECURSIVE  SOLUTION  OF  THE  STATE  TRAJECTORIES 

We  have  found,  in  Eq  .  (23),  the  solution  for  a  given 

perturbation  expansion  X.(t).  In  order  to  proceed  to  the 
next  order  (i+1),  we  neeJ  a  sampled  data  set  of  the  trajec¬ 
tories  of  X  (t).  This  will  permit  efficient  calculation 
of  the  Fourier  series  representation  of  the  inhomogeneous 
terms.  To  duplicate  the  procedure  outlined  above  for  each 
value  of  t  would  of  course  be  far  too  costly  in  terms  of 
time  and  computational  efficiency.  However,  an  alternative 
exists  as  a  result  of  the  particular  exponential  property 

A(k+l)At  AAt  AkAt  . 

e  =  e  e  (29) 


By  making  extensive  use  of  the  property  given  by  Eq  .  (29), 

we  can  produce  a  recursive  solution  of  Eq .  (23)  in  which 

each  data  point  (for  a  small,  fixed  time  interval  At)  is 
found  by  simple  matrix  multiplication  and  addition. 

It  can  be  shown  that  applying  the  property  in  Eq  .  (29) 

to  the  definitions  in  Eq .  (21b)  gives 


4j [ (k+l)At]  =  «  (At)$j(kAt)  (30a) 

*2it  (k+1)  At  ]  =  4  L  (At)*  (kAt)  +  *  (At)  4>3(kAt)  (30b) 


$-[ (k+1) At ]  =  4_(At) J  (kAt) 
3  3  3 


(30c) 


where  4j(0)=I,  an^  ^3(0), =1.  Now  let  us  define 

the  three  vectors  v,  .  to  be 

-J  1 


v j i [ ( k+ 1 ) At  ] 
v2  .  [  ( k+ 1 ) A  t  ] 
v3 . [ ( k+ 1 ) At  ] 


♦j [ (k+l)At]X.(0) 
42 . [ (k+1) At ]£0 
t  ( k+ 1 ) At ] £  q 


(31a) 

(31b) 

(31c) 


Evaluating  Eq.  (23)  at  t=[(k+l)At],  substituting  Eq.  (31) 
into  the  result,  and  by  substituting  Eq.  (31)  into  Eq.  (30) 
we  get 
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Xi[(k+l)At]  =  vli[(k+l)Atl  +  v  2  i  (  ( k+ 1 ) A  t  ] 


(32a) 
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—  i  i  (  (  k+  1 )  A  t  ]  =  <t1.(At)v1.(kAt)  ,  Vj.CO)  =  X^O)  (32b) 

—2  i [ ( k+ 1 ) At ]  -  ♦ j (At)v2i(kAt) 

(32c) 

+  $2 . (At)v3.(kAt)  X2i(0)  =  - 

v3i  [  (k+1)  At  ]  -  't3(At)  v3i  (kAt)  ,  3  ^  ( 0)  °‘Sl0  (32d) 

We  can  see  by  Eq .  (32)  that  we  need  only  to  calculate 
the  matrix  exponential  given  by  Eq.  (21)  once  in  order  to 
implement  this  recursive  solution.  Clearly  then,  the  solu¬ 
tion  of  each  order  requires  the  calculation  of  the  matrix 
exponential  twice;  once  at  t  =  tf  to  solve  for  the  initial 
costates  and  once  at  t*=At  to  set  up  the  recursive  procedure. 
When  the  solution  of  each  order  is  completed,  the  variables 
are  combined  according  to  Eq .  (8)  to  produce  the  open  loop 

optimal  trajectories.  We  again  stress  the  fact  that  the 
solution  to  the  nonlinear  optimal  control  problem  has  been 
produced  by  solving  a  sequence  of  strictly  linear,  constant 
coefficient  subproblems  without  resorting  to  any  iterative 
techniques.  We  shall  illustrate  the  effectiveness  of  the 
perturbation  method  with  numerical  examples  of  a  multi- 
axial  attitude  control  of  a  rigid  spacecraft. 

SPACECRAFT  EQUATIONS  OF  MOTION 

The  spacecraft  configuration  for  a  multi-axial  rigid 
body  is  immaterial,  but  we  ghall  assume  that  we  have  selected 
a  body  fixed  set  of  axes,  {b^},  located  a  the  mass  center 
of  the  body  and  aligned  with  the  principal  axes  of  the  space¬ 
craft.  Describing  the  orientation  of  the  body  fixed  frame 
relative  to  the  inertial  frame,  {  N^ .  }  ,  by  a  set  of  1-2-3 
Euler  angles  ( $ ,  0 ,  i|0  ,  we  can  show  that  the  angular  velocity 
of  the  spacecraft  in  body  fixed  axes  is  given  by 

w_  =  ul  — 1  +  u2  —2  +  w3— 3  (33) 

where 

•  » 

Wj  =  (+>c9cip  +  6sip  c(  )  =  cos(  ) 

u2  ■=  Gcip  -  $  c  0  s  ip  s(  )  =  sin(  ) 

•  • 

w  3  =  s  0  +  ip 
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Inverting  these  kinematic  equations,  we  obtain  the  Euler 
angle  rates  as 


• 

$  - 

U  C*/c0  -  IDjSlii/ce 

(34a) 

e  - 

(i)  j  s(Ji  + 

(34b) 

ij,  - 

-<jjjC^s9/c9  +  D2S<f)s9/c9  + 

(34c) 

As  indicated  previously,  we  must  expand  the  trigonometric 
functions  in  a  Taylor's  series.  In  doing  so,  we  shall 
neglect  all  quartic  and  higher  powers  of  the  variables, 
resulting  in  the  approximate  nonlinear  equations 


CM 

3 

1 

/“N 

CM 

CM 

CN 

CM 

d> 

s-' 

r—* 

3 

+ 

(35) 

2 

+  tdjlll  -  /2 

(36) 

-  w  ^  9  + 

(37) 

We  can  show  that  these  equations  are  adequate  to  describe 
maneuvers  of  j<20°  in  any/all  axes.  The  spacecraft  is  to 
be  controlled  by  three  external  control  torques  (u.  ,  i=l  ,2,3) 
applied  about  the  respective  principal  axes.  Therefore 
the  dynamics  of  the  multi-axial  motion  is  defined  by  Euler's 
equations  as 


U  1  7  1 1 

+ 

(I2 

_I3^w2U)37ll 

(38) 

U  2  7  1  2 

+ 

(I3 

(39) 

u  3  7  1  3 

+ 

(I1 

(40) 

where  I.  is  the  ith  principal  inertia.  The  motion  of  the 
spacecraft  is  governed  by  Eqs.  (35-40)  and  our  task  is  to 
find  the  control  torques  required  to  drive  the  spacecraft 
to  a  final  state  of  zero  attitude  and  zero  angular  velocity 
from  reasonable  arbitrary  initial  conditions  (defining  the 
fixed  target  orientation  to  be  the  inertial  frame). 


VV.V 


In  the  optimal  control  formulation,  we  shall  include 
the  control  torques  and  their  first  derivatives  in  the 
state  vector.  By  doing  so,  we  define  our  pseudo-control 
0.  =  ij  and  as  such,  we  shall  be  penalizing,  in  the  perform¬ 
ance  functional,  the  integral  norm  of  the  "control  accel¬ 
erations".  The  resulting  control  profiles  will  then  have 
prescribed  (usually  zero)  magnitude  and  slope  at  the 
initiation  and  completion  of  the  maneuver.  Therefore,  the 
augmented  state  and  control  vectors  are  defined  to  be 

£  [  $  6  ojj  u2  Uj  u 2  u3  Uj  u 2  u.j]T  (41) 


U  »  (Uj  u ^  u3] 


We  can  then  put  Eqs.  (35-40)  into  the  matrix  form  indicated 
by  Eq.  (1),  and  follow  the  procedure  outlined  in  the  develop¬ 
ment  of  the  perturbation  method. 

NUMERICAL  EXAMPLES 

We  demonstrate  how  well  the  perturbation  method 
performed  in  the  following  numerical  examples  by  summarizing 
results  from  integrations  of  the  exact  equations  of  motion 
(Eq.  (34)  and  Eqs.  (38-40)  ),  using  the  optimal  controls 
generated  by  the  perturbation  method,  and  comparing  the 
final  conditions  of  the  integrated  equations  with  the  desired 
final  conditions.  Repeating  this  procedure  for  solutions  of 
different  orders  in  the  expansion  given  by  Eq-.  (8)  demon¬ 
strates  how  the  accuracy  is  improved  with  each  order. 

CASE  1 

For  the  first  example,  a  small  angle,  rest-to-rest 
maneuver  in  2  seconds  is  performed  with  the  perturbation 
expansion  carried  to  second  order.  The  initial  conditions 


0.1  rad 


9  =  0.05  rad 


co  j  =  0 


=  0.15  rad 


Cl'  j  =  0 


We  select  the  weight  matrices  Q=0,  R  =  I  ,  and  S  =  1  .  F.  2  0  x  [1], 
The  large  value  of  the  matrix  S  serves  as  a  requirement 
that  the  final  conditions  be  rigidly  enforced.  In  this 
problem,  we  are  penalizing  only  the  "control  accelerations" 


and  the  final  state  errors.  Since  this  is  a  small  angle 
maneuver,  we  would  expect  the  solution  to  be  accurate,  since 
the  Taylor's  series  approximation  made  in  deriving  Eqs.  (35-37) 
will  be  very  accurate.  The  results  of  the  integration  of  the 
equations  of  motion  are  shown  in  Table  1  for  each  order.  We 
can  see  that  the  zero  order  solution  alone  is  quite  precise, 
as  we  suspected,  and  that  by  the  second  order  solution,  we 
have  obtained  nearly  two  orders  of  magnitude  improvement 
toward  the  final  conditions.  This'level  of  precision  is 
probably  all  one  would  require  for  open  loop  maneuver  controls, 
since  (i)  modeling  errors  and  implementation  errors  are  always 
present,  (ii)  if  high  precision  is  required,  we  would  always 
employ  a  terminal  feedback  control  for  precise  target  acqui¬ 
sition.  Clearly,  the  perturbation  method  is  providing  an 
acceptable  degree  of  accuracy  for  this  open  loop  problem. 


TABLE  1  FINAL  STATE  ERRORS 


Order  Approximation 


$ 

0 

-1  .  1  IE-3 

1 

-8 . 97E-5 

2 

-3 . 1 7E-4 

e 

3 . 53E-3 

2.87E-4 

-7 . 1  IE-6 

1  .  54E-3 

-5. 72E-5 

-1 . 08E- 5 

<*>i 

-4 . 48E-3 

7 . 22E-4 

.4 . 30E-5 

u:  2 

1  .  1  2E-2 

5 . 78E-4 

-4 . 48E-5 

w  3 

-  5 . 6  1  E  -  4 

7  .  72E-5 

9. 77E-6 

CASE  2 


In  this  example  we  shall  challenge  the  perturbation 
method  by  increasing  the  initial  attitude  to  larger  angles 
and  by  letting  the  spacecraft  be  tumbling  initially.  As 
in  the  previous  example,  we  seek  the  optimal  controls  that 
will  drive  the  system  to  a  zero  attitude  at  rest  with  a 
maneuver  time  of  2  sec.  The  initial  conditions  for  this 
case  are 

$  c  0.25  rad  w ^  =  0.4  rad/sec 

6  “  0.10  rad  Ui  2  =  0.2  rad  / sec 


= 


=  0.40  rad 


1.0  rad/sec 


The  weight  matrices  are  given  the  same  values  as  in  Case  1. 

In  like  manner,  the  results  of  the  perturbation  solution  to 
order  2  are  summarized  in  Table  2.  We  note  immediately  that 
the  zero  order  solution  is  ineffective  at  performing  the 
detumbling  maneuver  and  may  actually  aggravate  the  motion  of 
the  spacecraft.  However,  the  first  and  second  order  approxi¬ 
mations  converge  quickly  to  produce  acceptable  final  state 
errors  and  demonstrate  the  effectiveness  of  the  perturbation 
method.  We  would  conclude  that  the  control  trajectories 
produced  by  the  first  or  second  order  approximations  are 
indeed  adequate  approximations  of  the  open  loop  optimal 
controls  which  we  wished  to  determine.  The  trajectories 
of  the  Euler  angles,  angular  velocity,  and  controls  for 
the  second  order  solution  are  shown  in  Figs.  (1-3)  respec¬ 
tively. 


TABLE  2  FINAL  STATE  ERRORS 


Order  Approximation 


0 

1 

2 

4> 

-0.218 

-0.430E-1 

-0. 184E-1 

e 

0.315 

-0. 1 38E-1 

-0.435E-3 

U> 

0 . 679E-1 

-0. 123E-2 

0. 142E-2 

u  1 

-0 . 104E-1 

0. 157E-1 

0 . 740E-2 

OJ  2 

0.434 

0.2I8E-1 

-0. 152E-1 

0)3 

0 . 224E-2 

0  .  1  47E-2 

0. 161E-2 

CONCLUSIONS 

A  procedure  for  solving  nonlinear,  open  loop,  optimal 
control  problems  has  been  presented.  In  this  approach, 
an  asymptotic  perturbation  method  is  applied  to  the  problem 
thereby  eliminating  the  traditional  dependence  on  iterative 
numerical  methods.  The  nonlinear  system  is  "separated" 
into  a  set  of  nonhomogeneou s ,  linear,  optimal  control  prob¬ 
lems  that  may  be  solved  sequentially  as  "independent" 
systems.  Upon  combining  the  solutions  of  the  subproblems 
in  a  straightforward  power  series,  an  optimal  control  for 
the  nonlinear  system  is  generated.  This  novel  process  for 
solving  nonlinear  optimal  control  problems  is  a  result  of 
the  marriage  of  a  simple  analytical  technique  (the  pertur¬ 
bation  method)  and  a-apowerf  ul  numerical  algorithm  (the 
matrix  exponential). 


We  have  applied  the  perturbation  method  to  several 
classes  of  problems  (Ref.  7,8)  and  have  found  it  to  be 
most  effective  for  a  large  family  of  nonlinear  problems, 
including  large  degree  of  freedom  systems.  The  material 
presented  in  this  paper  is  a  further  evaluation  of  the 
method  in  its  present  format  regarding  nonlinear  three 
dimension  spacecraft  control.  We  are  continuing  to  improve 
the  method  in  an  effort  to  achieve  greater  numerical  effi¬ 
ciency  while  simultaneously  applying  the  approach  to  a 
greater  number  and  wider  variety  of  nonlinear  problems. 
Higher  order  solutions  to  common  problems  can  be  completed 
quickly  and  with  minimal  programming  effort  in  a  manner 
that  cannot  be  matched  by  purely  analytical  methods,  yet 
with  an  accuracy  that  is  consistent  with  other  numerical 
methods.  We  anticipate  that  the  quasi-analy tical  pertur¬ 
bation  method  will  prove  the  most  useful  approach  for 
solving  many  nonlinear  problems. 
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plus  boundary  conditions:  zxo)  =  zn,  p(tf)  =  Sz(tf)  or  Z(tf) 


PRELIMINARIES 


Two  Point  Boundary  Value  Problem 


Due  to  the  nonlinear  terms  in  h(Xjt),  exact  analytical  solutions  ofEq .  (4)  are 
most  often  impossible .  A  variety  of  iterative  techniques  are  available;  they 
are  often  expensive  due  to  initial  ignorance  of  a  " good  starting  estimate  (of 
P(t0 )  or  p(tj) )  required  for  reliable  convergence.  We  seek  to  avoid  iteration 

through  use  of  a  perturbation  method  &  " quasi-analytical 99  integration.  »> 


But...  how  do  we  make  efficient  algorithms?  Does  convergence  occur  in  the  " real  world"?  Can  the 
above  be  implemented  in  a  way  which  automates  the  algebra  usually  associated  with  perturbation 
methods?  What  about  secular  terms?  Does  this  approach  apply  to  systems  of  non-trivlal  dimen¬ 
sions  &  "messy"  nonlinear  terms?  We  have  made  some  progress  In  answering  these  questions. 
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Subject  to:  z  =  Az  +  Du.  The  necessary  conditions  have  the  identical 
form  as  those  developed  in  the  foregoing.  Penalizing  the  control  derivatives 
has  been  found  most  constructive  in  frequency-shaping  the  torque 
profiles  to  decrease  excitation  of  the  poorly  modeled  higher frequency  modes. 
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Case  1  Optimal  Detumble /Attitude  Acquisition  Maneuver 


Case  2  Numerical  Results  for  the  TPBVP  Solution 


(out-of -plane)  0.317E-2  0.29 
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K  g 

s  « 

■g  8  I 

«  ■«  2 

Js  to  *5: 

S  S3  g 

«  SS, 
T*  *■  ■s* 

Sf  «  h 

K  a  a 

a  .5 

S  ©  *3 

"2  *  § 

s  8  £ 

o  *  a* 

^  8  8 

^  a 

•*»  *5; 

a  a 

I  g 

a,  ft 

c>  c 

ft  fto 

Sk. 

ft  Q 
ft 


ft  Jo 

*  3 

.1  g 

|a 

k  't* 

f!  g 

&,  §, 
6  8 

§  •« 
•a  -h 
£  ^ 
g  s 

«  g 
s  ft 

* 

S  ^ 

£  « 

ft  ft 
ft  ft 
ft  nJ 

<■>  t: 
*8  § 

^  •a 

5  5 


.9  A 

4  | 

*  g 

•ft  *f* 

tO  k 

a  S> 

ft, 

£T^ 

ft  ft 

2 
k  J* 

«  % 
ft  *ft 

ft  ^ 

a  \'  ft 

I  fc’g 

r|  I 

ft  ^  -a 

ft  ^  ^ 

to  a  g, 

a  vs  -2 

a  a  ft 
a,  ft  >• 

gf  ft  .ft 
^  a  « 
J*  ft,  60 
•a  h*  ft 
.ft  ft  a 

a -8 1 

ft  a  ft 

8  2  E 


•a  "a 
•1  ^  ^ 
*  £  g 
5  *1  a 

•a  5  5 

•a  *o  ft 

■*»*  ft  *ft 

'S' -I  § 

|*s 

a  lit 

a  o,  k 

ft  C  ft 

a  ft  *a 

ft  p>  k 

a  a  a 

k  ft  ft*. 

ft,  a  a* 
ft  ft  so 

«  s  g 

C  •S  3 

6,1  & 

ft,  V  •** 

8  &  8 

g  •$  1 

o  a  a 
a  -ft  ft 
a  ^  ^ 

ft  « 

vs  ,&>  >, 
ft,  C  "a^ 
g*  k  a 
a  «  «a 

^  ft  ft 

ft^  a  vs 

>,  as  !§ 
ts  ft  h 

•2:  ft  ft 
a 


fc  g 

ft,  *a 

%>  * 

S  *8 
'S'l 

ft  ft 
ft  ft 

5  § 

£?  2 
ft  ft 

p*  vs 

I  s 

fe  Si 

1  =3“ 

*  ‘B  « 

2  ft  vs 

ft  *C  ft 
■8  1  1 
s'  i  I 

a  k  ft 
c  a  ft 

.2  a  *2 

3  .2  ■$ 
o2^ 
ft  a  k 

ft  a  ft 

•2  §”S 
3  *  tc 

^  ft  ft 

M  a 

111 

S  3  2 


s  g 
•§ * 
8  ■§ 

■8  * 
Si  "2j 

«  .a 


ft  .t 

3  .S 

8  g 

•s  ?. 

ft  ft 
p.  ft 

•  *  *m 

•«*  -v. 

ft  a 


a  a 

•a  ? 

■8  | 
VS  to 
ft  k 

5  <* 

ft 


AD-A172  716  OPTIMIZATION  OP  CLOSED  LOOP  EIGENVALUES  MANEUVERING  2/2 
VIBRATION  CONTROL  AN  <U>  VIRGINIA  POLVTECHNIC  INST  AND 
STATE  UN IV  BLACKSBURG  DEPT  OF  E  J  L  JUNKINS 
UNCLASSIFIED  11  MAV  8b  AFOSR-TR-86-0905  F49620-81-K-0012  F/G  22/2  NL 


Attachment  #6 


<$ 


i5  Identification  of  Vibrating  Flexible  Structures 

_  S.  Rajaram  and  J.L.  Junkins 


I 


LVV 


VOL.  8,  NO.  4,  JULY-AUGUST  1985 


J.  GUIDANCE 


463 


Identification  of  Vibrating  Flexible  Structures 

S.  Rajaram*  and  J.L.  Junkins+ 

Virginia  Polytechnic  Institute  and  Slate  University,  Blacksburg,  Virginia 

This  paper  presents  novel  identification  schemes  to  determine  model  parameters  of  vibrating  structures.  A 
time-domain  identification  method  using  transient  response  is  discussed  first.  Next,  a  steadi-stale  response 
method  using  nonresonant  harmonic  excitations  is  considered.  An  especially  attractive  method  for  unique's 
identifying  the  parameters  of  a  structure  using  both  free  and  forced  response  is  also  discussed.  Numerical  results 
show  that  the  methods  are  relatively  immune  to  the  presence  of  damping  and  mans  loss-frequent  *  modes  with 
repealed  or  closely  spaced  frequencies. 


Introduction 

CT1VE  control  of  large  space  structures  (l.SS) 
necessitates  a  sufficiently  accurate  estimate  of  the 
parameters  so  that  control  laws  can  be  tuned  on-orbit  to  en¬ 
sure  stability  and  permit  l^s  control  effort  to  be  expended. 
Algorithms  for  design  of  insensitive  or  adaptive  controls  are 
not  attractive  due  to  the  large  number  of  degrees  of  freedom 
to  be  controlled.  In  fact,  for  most  LSS  application,  the  only 
feasible  approach  appears  to  be:  I)  identify  the  structural 
parameters  then  2)  use  this  information  to  adjust  the  gains, 
and  perhaps  3)  use  adaptive  methods  to  change  a  small 
number  of  critical  parameters  in  real  time.  This  paper  ad¬ 
dresses  issue  1  above. 

Transient  Response  Identification  Method 

Many  structural  modal  identification  methods  are  available 
which  extract  modal  characteristics,  i.e.,  natural  frequencies 
and  mode  shapes  from  a  set  of  resonant  steady-state  responses 
due  to  a  large  number  of  harmonic  excitations.  These  methods 
encounter  analytical  and  numerical  difficulties  when  the 
system  frequencies  are  closely  spaced  and  the  “single  mode 
resonant  response"  assumption  is  used.  Also,  the  time  re¬ 
quired  to  achieve  steady  state  may  be  prohibitively  long  for 
lightly  damped,  low-frequency  structures.  Time-domain 
techniques' 4  for  structural  identification  were  first  proposed 
by  Ibrahim.  Ibrahim's  time-domain  (1TD)  method  is  a  modal 
identification  scheme.  The  ITD  method  has  been  successfully 
applied  to  reduce  measurements  from  several  laboratory  ex¬ 
periments,  however,  this  method  has  been  found  to  lack 
reliable  robustness.  In  some  applications,  rank  deficient  linear 
systems  are  encountered.  Recently,  Juang  and  Papa'  have 
developed  a  more  robust  time-domain  modal  identification 
method;  this  method  is  based  upon  judicious  use  of  singular 
value  decomposition. 

Identification  in  Configuration  Space 

Consider  a  vibrating  structure  governed  by  the  following 
linear  matrix  differential  equation: 

Mx  +  Cx  +  Kx  =f  (I) 

where  x  is  the  n  x  I  configuration  vector  of  physical  displace¬ 
ment,  /  the  n  x  I  force  vector,  M  the  nxn  symmetric  positive 
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definite  mass  matrix,  C  the  nxn  symmetric  positive 
semidefinite  damping  matrix,  and  A'  the  n  x  n  symmetric 
positive  semidefinite  stiffness  matrix.  Dots  denote  differentia¬ 
tion  with  respect  to  time. 

It  is  assumed  that  a  LSS  can  be  satisfactorily  modeled  in  (he 
form  given  by  Eq.  (1).  Our  objective  herein  is  to  identify  the 
poorly  known  coefficient  matrices  M,  C.  and  A  or  some 
parameterication  therefore,  e.g.,  the  system  eigenvalues  and 
eigenvectors.  Equation  (I)  can  be  rewritten  as 

M 

|i'(f)  x'(f)  jt'  (/)  C  =  /'(()  (2) 

:  A 

T  denotes  the  matrix  transpose  operation. 

Now,  consider  an  idealized  measurement  process  wherein 
the  position,  velocity  acceleration,  and  forces  arc  measured  at 
discrete  instants,  sav  i,,l, . Upon  writing  m  measure¬ 

ment  equations  (»/>3n)  identical  to  Eq.  (2),  one  for  each 
measurement  time,  the  resulting  matrix  equations  can  be  writ¬ 
ten  as 

z.p  =  u  (3i 

where  Z  is  an  m  x  3  n  coefficient  matrix,  whose  j  row  contains 
measurements  of  the  system  response  at  lime  I/. 

yth  row  of  /  -  |i'  (/,  )x'  < /,  )Jr'  ( t, )  |  (4) 

G'is  a  m  *  n  matrix  containing  the  following  forcing  functions: 

yth  row  of  U  =/'(/, !  (5) 

P-  |.V  (  A  | '  (6| 

P  is  a  3n  x  n  matrix  containing  the  unknown  mass,  damping, 
and  stiffness  parameters. 

Since  the  number  of  elements  in  each  column  of  P  is  Jn  and 
m>3n,  Eq.  (3)  ovcrdctermincs  the  columns  of  P.  The  least 
squares  solution  for  P  is  given  as 

P^LU  <7| 

where  the  least-squares  operator  is  formally 

/  =  (/'/)  (X) 

for  large  systems,  of  course,  the  explicit  inversion  should  be 
avoided  in  favor  of  the  Q  P  reduction,  C'holesky  decomposi¬ 
tion,  or  a  singular  value  decomposition  approach  lor  in 
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creased  efficiency  and  robustness.1’  The  computations  subse¬ 
quently  summarized  in  this  paper  were  done  using  the  Q-R 
algorithm.  The  only  theoretical  requirement  is  that  the  least- 
squares  coefficient  matrix  Z  hare  full  rank  ( 3n ).  Physically, 
this  rank  condition  can  be  achieved  only  if  all  degrees  of 
freedom  participate  in  the  response.  Hence,  a  fundamental  re¬ 
quirement  for  identifying  the  structure  is  that  the  excitation 
should  hare  sufficient  energy  and  frequency  content  to  excite 
the  higher  modes  of  the  system.  Qualitatively,  it  is  also  evident 
that  the  actuator  locations  and  phase  distribution  of  the  ac¬ 
tuator  input  are  likely  to  be  important.  Thus,  the  mass,  damp¬ 
ing,  and  stiffness  matrices  of  the  structure  can  be  identified,  at 
least  in  principle.  The  method  is  obviously  straightforward. 
However,  it  requires  that  the  number  of  forces  equal  the  order 
of  the  system.  Also,  acceleration,  velocity,  and  displacement 
are  to  be  measured  at  all  of  the  degrees  of  freedom.  These  re¬ 
quirements  pose  obvious  practical  difficulties.  It  is  shown  in 
Ref.  7  that,  in  order  to  use  a  smaller  number  of  forces  than  the 
number  of  degrees  of  freedom  of  the  system,  a  priori 
knowledge  of  the  mass  matriv  is  required.  It  should  be  noted 
that  the  method,  as  presented  above,  does  not  require  or  ex¬ 
ploit  the  symmetry  of  the  system  matrices  and,  hence,  is  ap¬ 
plicable  to  a  general  dynamic  system  involving  gyroscopic  and 
circulatory  forces.  It  is  also  evident  that  the  size  of  the  linear 
systems  which  must  be  solved  is  in.  Therefore,  unless  the 
matrices  do,  in  fact,  possess  special  properties,  it  is  anticipated 
that  practical  computational  restrictions  will  require  n<  50  for 
this  approach.  Of  course,  the  heavy  redundancy  implicit  in  M, 
C.  and  A  as  descriptions  of  distributed  mass,  damping,  and 
stiffness  often  can  be  eliminated  in  terms  of  fewer  physical 
parameters,  but  the  estimation  process  then  must  be  coupled 
with  the  structural  modeling  (e.g.,  finite  element)  process. 


The  matrices  M  'A',  M  'C.  and  B,  can  be  determined, 
following  a  procedure  analogous  to  that  outlined  previously 
for  configuration  space  identification.  It  is  evident  from  Lq. 
(14)  that  the  least-squares  coefficient  matrix  includes  the  force 
vector.  Immediately,  it  can  be  inferred  that  the  force  rector 
should  form  an  independent  set  for  unique  identification  of 
the  system  parameters.  This  statement  holds  true  for  con¬ 
figuration  space  identification  also.  It  can  be  observed  that  jr 
in  the  least-squares  coefficient  matrix,  viz.,  Eq.  (3),  becomes 
independent  of  other  variables  only  through  the  forcing 
functions. 

Identification  Using  Orthogonal  Polynominals 

The  measurement  of  acceleration,  velocity,  and  displace¬ 
ment  at  all  degrees  of  freedom,  as  required  by  the  methods 
presented  earlier,  poses  obvious  practical  difficulties.  To  ob¬ 
tain  partial  relief  from  this  requirement,  an  orthogonal  iden¬ 
tification  scheme  is  proposed.  Orthogonal  polynomials  can  be 
used  to  represent  completely  any  function  to  a  required  degree 
of  accuracy.* 

Consider  the  lower  portion  of  the  state  equations,  viz.,  Lq. 
(13).  Also,  it  is  assumed  that  the  accelerations  and  forces  are 
measured  at  discrete  instants  of  time.  Then  they  can  be  ex¬ 
panded  using  orthogonal  polynomials  such  as  Chebvshev, 
Legendre,  etc. 

x  =  P,T(t)  (15) 

where  P ,  is  a  rectangular  coefficient  matrix  and 

T(t)  =  |  T0U)T,U)..  7\_  / (/)  1 r  (16) 


Identification  in  State  Space 

For  control  applications,  the  system  dynamics  is  expressed 
in  state  equations.  Introducing  the  ”2n”-dimensiona)  state 
vector 

*</)=  [*'(/)  *'(/)]7  (9) 

Equation  (1)  can  be  written  as 


is  a  column  vector  consisting  of  orthogonal  polynomials.  The 
intergral  of  T,  (t)  can  be  expressed  via  a  recurrence  relation¬ 
ship  involving  T, and  T,  /.  Integrating  Eq.  (15) 

*  =  P;T(l)  +  C'  =)'i  4  C'  (17) 

Integrating  further 

jr=  P,TU)  +  c't  +  c"  =y+  c'ts  c"  (IS) 


*  -  rig  v  Bf 


M  A  -  M  ( 


is  the  plant  dynamic  matrix  and 


where  P,.  P~,  and  P,  arc  n  x  \'  matrices  containing  the  expan¬ 
sion  coefficients.  By  substituting  Eqs.  (15).  (17),  and  (IS)  into 
Eq.  (13),  the  least-squares  problem  can  be  constructed  as 

i  ~(M  1 K)'  } 


l  .»•'(/)  r'it)  fit)  l  1 1 


(A/  'O' 


jf'(/>  (14) 


is  the  control  distribution  matrix.  The  structure  o!  H  is  depen¬ 
dent  upon  type  and  location  ol  the  force  inputs.  When  all 
degrees  of  freedom  are  not  excited,  the  force  vector  contains 
zero  entries  and  H  will  be  a  (2  n>n.)  reclangulai  matrix, 
where  n ,  is  the  number  of  excitations.  The  unknown 
parameters  to  be  identified  are  the  elements  ol  matrices  f  and 
B 

Consider  the  lower  partition  of  l  q  (10) 

x  (/)  -  M  'hxtn-M  '(’xtt)  •  />,/(  / )  (13) 


I  quation  (13)  can  be  written  as 


(  -  M  A  ) '  1 


1  tn  |jr'(/l  Jr'lri  /'(/)!  (  -  '/  (  ) ' 


M  Ac 


A 1  'Ac'-Af  'Cc 


Af  'A.  Af  'C.  P,,  d,.  and  d.  can  be  estimated  Irorn  Eq.  ( 14) 
c  and  c"  can  be  determined  using  Eqs.  (20)  and  (211  Thus, 
the  numbet  ol  measurements  are  reduced  by  a  lactot  ol  •  and 
the  initial  displacement  and  velocity  vectors,  usually  the 
equilibrium  positions  ol  the  structure,  are  also  estimated  along 
with  the  parameters,  it  is  also  possible  to  use  velocity  oi 
displacement  measurements  alone  Although  several  oi 
thogonal  polynomials  exists,  the  use  ol  Chebvshev 
polynomials  have  found  a  wide  application'  in  solving  lineat 
and  nonlinear  diffcicnlial  equations  Since  we  are  concerned 
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with  the  inverse  problem  <i.e. ,  given  the  response,  the  best 
estimate  of  the  system’s  parameters  is  to  be  determined), 
Chebyshev  polynomials  seem  to  be  the  natural  choice. 

Numerical  Results  for  Identification 
from  Transient  Free  Response 

Four  specific  linear  systems  are  considered  as  representative 
examples;  a  spring-mass  damper  system  (Fig.  I),  a  plane  truss 
(Fig.  2),  a  cantilever  beam  (Fig.  3),  and  a  rectangular  mem¬ 
brane  (Fig.  4)  are  considered  to  study  the  effects  of  I )  repeated 
low  frequencies  and  rigid-body  modes,  2)  damping,  3)  choice 
of  excitation  and  number  of  excitations  vs  degrees  of  freedom, 
4)  excitation  frequency  vs  system  natural  frequency,  5) 
measurement  errors,  measurement  duration,  and  sampling  in¬ 
terval,  and  6)  model  truncation  errors.  Synthetic  measured 
data  are  generated  for  each  case  using  known  parameter 
values.  Table  1  gives  the  undamped  eigenvalues  of  the  ex¬ 
amples.  The  proposed  identification  schemes  performed  very 
well  for  all  four  examples.  The  results  are  summarized  below . 

The  plane  truss  example,  as  can  be  seen  from  Table  I,  has 
three  repeated  eigenvalues  and  three  zero  eigenvalues  cor¬ 
responding  to  the  rigid-body  modes.  Arbitrary  viscous  damp¬ 
ing  is  included  in  the  first  and  second  example  problems.  It  is 
found  that  the  identification  algorithms  recovered  the  system 
matrices  without  any  difficulty.  It  is  clear  that  restrictive 
assumptions  such  as  proportional  or  negligible  damping  are 


unnecessary  for  the  success  of  the  identification  algorithms. 
Three  excitation  types  were  considered:  1 )  harmonic,  2)  bang- 
bang  (rectangular  wave),  and  3)  frequency  swept  harmonic 
(harmonic  excitations  with  time-varying  frequency).  Har¬ 
monic  excitations  yield  good  results  for  the  spring-mass 
damper  and  cantilever  beam.  For  the  plane  truss,  bang-bang 
and  frequency  swept  harmonic  excitations  are  useful.  Or¬ 
thogonal  polynomial  identification  of  the  plane  truss  with 
bang-bang  excitation  did  not  recover  the  parameters  very  well. 
The  orthogonal  polynomials  are  unable  to  represent  the  ac¬ 
celerations  satisfactorily.  However,  with  frequency  swept  har¬ 
monic  excitation,  the  orthogonal  identifier  recovered  all 
parameters  accurately.  Thus,  the  orthogonal  identification 
scheme  apparently  works  best  w  ith  smooth  and  continuous  ex¬ 
citations.  The  number  of  excitations  required  for  the  iden¬ 
tification  of  spring-mass  damper  and  cantilever  beam  can  be 
as  few  as  one.  For  the  plane  truss,  a  minimum  of  five  excita¬ 
tions  are  needed. 

Upon  varying  the  frequency  of  excitation  over  the  range  of 
the  system  natural  frequencies,  no  significant  degradation  in 
the  performance  of  the  algorithm  is  observed.  Measurement 
errors  introduce  estimate  errors,  of  course.  The  effect  of 
discrepancies  between  commanded  and  realized  excitations  is 
also  studied  by  including  noise  in  the  excitations. 

The  cantilever  beam  and  a  simply  supported  membrane  are 
chosen  to  illustrate  the  effect  of  model  truncation.  These  are 
distributed  systems;  we  are  interested  in  obtaining  a  discrete 
representation.  The  response  of  the  cantilever  beam  is 
obtained  by  using  the  eigenfunctions  and  assuming  that  six 
modes  participate  in  the  response.  Hence,  a  sixth-order  model 
is  obtained  first.  The  model  identified  using  either  a  bang- 
bang  or  harmonic  excitation  is  the  same.  A  reduced-order 
model  (fourth  order)  is  identified  next.  Table  2  compares  the 
exact  eigenvalues  with  those  obtained  from  identified  models. 
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Fig.  3  Cantilever  beam,  mix)  is  the  mass  unit  length.  2.4  kg  m:  / 7t.il 
the  bending  stiffness,  495  \/m:;  /  the  length.  3.0  m;  •  the  position  of 
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Fig.  I  Spring-mass  damper. 
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Fig.  2  Plane  truss. 
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Fig.  4  Rectangular  membrane.  *  is  the  position  of  sensors,  and  • 
the  position  of  actuators. 


Table  I  Eigenvalues  of  example  problems 


Spring  mass 

Plane 
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ll  is  also  noted  that  the  harmonic  excitation  resulted  in  a 
model  that  fits  the  measured  acceleration  fairly  well.  This  is 
significant  since  accelerometers  are  the  most  commonly  used 
sensors  for  vibration  measurement.  The  bang-bang  excitation, 
being  rich  in  harmonic  content,  is  able  to  excite  the  higher 
modes  considerably  and  thus  affects  the  identification  of  the 
reduced-order  model.  Table  3  gives  the  results  for  the 
membrane. 


Identification  Using  Nonresonanl 
Harmonic  Excitations 

The  use  of  harmonic  excitations  for  the  identification  of 
vibrating  structures  have  received  the  attention  of  several 
investigators. * 14  Raney1'  used  such  a  scheme  to  successfully 
identify  the  effective  masses,  stiffnesses,  and  damping  for  a 
lightly  damped  structure  having  widely  separated  modes;  the 
structures  studied  were  I  10-  and  I /40-scale  models  of  the 
Saturn  launch  vehicle.  Several  methods  are  suggested1' 14  to  ex¬ 
tract  the  normal  modes  from  measured  response.  However, 
the  use  of  normal  modes  is  questionable  when  the  damping  in 
the  system  is  not  a  proportional  type.  The  methods  using  reso¬ 
nant  harmonic  excitations,  as  mentioned  earlier,  encounter 
both  experimental  and  computational  problems,  when  the 
system  frequencies  are  closely  spaced.  A  novel  identification 
scheme  using  nonresonant  harmonic  excitations  is  presented  in 
this  section.  The  proposed  scheme  differs  significantly  from 
several  of  the  existing  methods  which  generally  use  resonant 
harmonic  response  to  obtain  the  model  parameters.  The 
method  requires  that  the  structure  be  damped;  the  damping, 
however,  can  be  arbitrary  viscous  damping. 


Configuration  Space  Identification 

Once  again  consider  Eq.  (I ).  l  et  the  excitation/)/)  be  given 
as 


u,k sinful/  +  ou  ) 
a,jsin(u;t/  +  o  -i ) 


/.I/)-- 


a,.i  sin  (,*•,/-*  o...  ) 


Table  2  Comparison  of  eigenvalues  for  cantilever  beam 


1  \»Il  t 

I  igctnaiucs  obtained  f  rom 

Sixth-order 

model 

I’Oiitih'ordcr  model 

Harmonic 

Bane  bane 

5.6105 

5.6545 

5.6282 

It  2121 

15.160" 

0" 

55.1664 

55.2092 

56.0591 

98  4M 

98.4524 

98.5540 

98.5429 

192.9246 

192  9292 

192.9285 

192.9575 

.MX  9IK2 

5  IK. 8*50 

4^6.40X6 

4"  5 .9985 

Table  3 

Comparison  of  eigen* allies  for  rectangular  membrane 

1  iccnvalucN  obtained  from 

1  oiittli  order 

model 

Ixaet 

model 

Harmonic 

Hang  bane 

21  0028 

21 .0008 

20.49 7  6 

21  22S4 

41  9l't 

41  9201 

42.070*7 

41  9604 

2i  pm 

21  0860 

21  42*6' 

51.0774 

42  00«« 

41  9602 

21  46  "6 

21  4405 

21  4256-' 

54.5880 

42.1521 

42  l  491 

■*  \ 

p i* r  t  |!tvcn  b\  Otitib 

. attic  totll l IK  IX  IH'j.'kxICll 

or 


fk  (!)  =.9*sinu)A/  +  Cjcosuij/  k  =  l,2 . m  (22) 

Sk  and  Ck  are  the  amplitudes  of  the  sine  and  cosine  com¬ 
ponents  of  the  excitation.  The  steady-state  response  of  the 
structure  then  can  be  written  as 

x(t)  =  Asmu\l  +  floosie*/  k  =  l,2 . m  (23) 

The  structure  is  subjected  to  “m"  harmonic  excitations  at  fre¬ 
quencies  ui,,us, . The  excitation  frequencies  can  be 

chosen  arbitrarily  and  need  not  coincide  with  the  system  fre¬ 
quencies.  For  each  excitation  frequency  “u >k,"  the  steady- 
state  amplitudes  Ak  and  Bk  of  the  displacement  are  measured. 
Using  Eqs.  (22)  and  (23)  in  Eq.  (2),  we  form  the  matrix 
equation 


ZP=U  (24) 

where  P  is  the  same  as  in  Eq.  (6). 

A'th  row  of  Z=  [ -wiAl  -ukBl  A[  j  (25) 

(A  +  m)th  row  of  Z  =  [ B\  J  (26) 

Cth  row  of  U-Sl  (27) 

(k  +  rw)ih  row  of  U=  Cl  (28) 


For  m> 3/1/2,  Eq.  (24)  represents  an  overdetermined  system 
of  equations.  The  least-squares  solution  for  P  is  given 
formally  as 


P=(Z'Z)  'Z’U  (29) 

Thus,  system  matrices  M,  C,  and  A'  can  be  identified  directly 
from  steady-state  response.  The  same  response  data  can  be 
used  to  identify  the  system  in  state  space.  Also,  the  amplitudes 
of  displacement  must  be  measured  at  every  degree  of  freedom. 
Alternatively,  amplitudes  of  accelerations  can  be  used.  In  that 
case,  Eqs.  (25)  and  (26)  become 

Ath  row  of  Z=  \A\  B’k/u:k  -.4, 7  Ac;]  (30) 

(*'  +  »<)th  row  of  Z=  \B[  -  A  [  /uik  (31) 


Ak  and  Bk  are  the  amplitudes  of  acceleration.  It  is  assumed  in 
the  subsequent  discussions  that  displacement  amplitudes  are 
measured  and.  hence,  consider  Eqs.  (25)  and  (26)  only. 


State-Space  Identification 

Consider  once  again  Eq.  (13),  which,  using  Eqs.  (22)  and 
(23),  becomes 


1  -w>47  j 

■l,7 

■v; 

(  -  A/  -A)'  ! 

1 

I  -4. -I7,  j  _ 

a;„  -u 

s!„ 

(  -M  ()' 

«'  -,.4' 

9 

Pi 

-  u-;. 

B'  u-,,.47, 

c!„ 

(32) 

Equation  (32)  can  be  solved  via 

least  squares  to  obtain  \f  'A. 

M  'C,  and  />,.  It  should  be  noted  that  the  amplitudes  of  the 
excitations  directly  enter  into  the  least-squares  coellicient 
matrix.  Hence,  they  should  form  an  independent  set  This  is 
achieved  by  varying  the  phase  ol  the  excitation,  i.c..  o,,, 
0.k . 0„,  in  Eq  (22)  in  a  nonlinear  fashion.  This  is  precisely 
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the  reason  for  choosing  the  excitation  of  the  form  given  by  Eq. 
(22).  The  same  requirement  holds  true  for  the  configuration 
space  identification,  since  the  amplitudes  of  Alk  and  Blk  can 
be  linearly  independent  only  via  the  excitation  amplitudes. 


Numerical  Results  for  a  Steady-Slate 
Response  Identification 

Two  numerical  examples,  viz.,  a  spring-mass  damper  (Fig. 
1)  and  a  plane  truss  (Fig.  2),  are  considered  to  study  the  effects 
of  several  implementation  issues  such  as  1)  closely  spaced/ 
repeated  frequencies  and  rigid-body  modes,  2)  number  of  ex¬ 
citations,  and  3)  choice  of  excitation  frequency.  The  results 
are  summarized  under  each  of  the  items  (1-3). 

The  plane  truss  example  is  used  to  study  the  effects  of 
repeated  frequencies  and  rigid-body  modes.  Thirty  excitation 
frequencies  were  used  ranging  from  0.6  to  0.89  rad/s  in  steps 
of  0. 1  rad/s.  All  of  the  parameters  are  identified  exactly  in  the 
absence  of  measurement  noise.  Thus,  the  proposed  scheme  is 
capable  of  identifying  systems  with  closely  spaced  frequencies. 
For  state-space  identification,  the  number  of  excitations  can 
be  as  low  as  three  for  the  spring-mass  damper  system. 
However,  the  plane  truss  requires  a  minimum  of  four 
excitations. 

The  excitation  frequency  was  varied  within  the  range  of  the 
system’s  undamped  natural  frequencies.  The  authors  con¬ 
sidered  those  cases  where  excitation  frequencies  were  below 
the  lowest  frequency  of  the  structure  and  above  the  highest 
frequency  of  the  structure.  No  significant  degradation  of  the 
identification  scheme  was  observed,  thus  validating  the  fact 
that  the  excitation  frequencies  can  be  chosen  fairly  arbitrarily. 


Identification  Using  Free-  and 
Forced-Response  Data 

Identification  of  mass  and  stiffness  matrices  from  model 
test  results  has  been  reported  by  several  authors.  The  objective 
of  these  identification  schemes  is  to  modify  an  a  priori  mass  or 
stiffness  matrix  so  that  measured  eigenvalues  and  eigenvectors 


agree  with  those  of  the  analytical  model.  Berman,' 


minimization  procedure,  developed  a  noniterative  scheme 
based  on  the  orthogonality  relationships  of  eigenvectors  for 
computing  a  "nearest  neighbor”  update  of  the  mass  matrix. 
Following  Berman's  approach,  Wei'”  developed  a  related 
method  for  correcting  the  stiffness  matrix.  Chen  and  Wada19 
discuss  an  interactive  system  parameter  refinement  procedure, 
employing  the  Jacobian  matrix  (consisting  of  the  derivatives 
of  eigenvalues  and  eigenvectors  with  respect  to  system 
parameters).  Recently,  Chen  et  al.-"  applied  a  first-order 
matrix  perturbation  approach  to  identify  the  mass  and  stiff¬ 
ness  matrices.  Other  related  approaches  can  be  found  in  Refs. 
21  and  22. 

Free-rcsponsc  data  are  used  herein  to  estimate  the  eigen¬ 
values  and  eigenvectors  of  the  system.  Using  orthogonality 
conditions,  the  matrices  equal  to  system  matrices  multiplied 
by  unknown  scale  factors  arc  determined  initially.  These  scale 
factors  are  then  uniquely  estimated  by  subjecting  the  system  to 
known  forces  and  measuring  the  acceleration,  velocity,  and 
displacement  at  several  locations.  The  approach  presented 
herein  embodies  a  fundamental  advantage:  perfect  measure¬ 
ments.  lead,  to  within  truncation  errors  and  arithmetic  errors, 
to  the  true  system  parameters. 


Berman's  Method:  A  Summary 

In  the  absence  of  damping  F.q  (I)  reduces  to 


Mx  r  A  jt  -  0 


The  orthgonahty  conditions  of  the  system  described  by  F.q. 
(33)  arc 

E1  ME  -  I  (34a) 


E'  KE  -  diag ( teywy . mi )  *  |  w’  I 


E:  (n  x  n)  modal  matrix 


/:  (n  xn)  identity  matrix 


w/.ui, . u!„  are  the  natural  frequencies.  Note  that  in  Lqs.  (34) 

the  eigenvectors  are  normalized  with  respect  to  the  mass 
matrix  so  that 


ejMe,  =  l  i=l,2 . n 


e,  is  the  ith  eigenvector  (/th  column  of  /-.').  It  will  be  assumed 
that  the  measured  modal  matrix  is  square. 

Berman"'  assumes  an  analytical  mass  matrix  M,  The 
measured  eigenvectors  are  normalized  with  respect  to  M ,  so 
that 


e]M.xe,  =  l  i=l,2 . n 


Letting  AM  be  the  desired  correction  matrix,  Berman 
minimizes  the  Euclidian  norm 


<  =  l N  '  AM  \  'I 


subject  to  the  constraint  equation 


E'ME=l-mu  M  =  M ,  -t  AM 


where  mu  =  E1  M .,£'  is  a  nondiagonal  matrix  having  unity  as 
diagonal  elements.  Choosing  N  =  M.t  as  the  weight  matrix. 
Berman  obtains 


AM  =  MxEma  1  (I-md)mu  'I  'M  , 


It  can  be  seen  from  Eq.  (38)  that  AM  is  symmetrical  and  deter¬ 
mined  to  satisfy  the  orthogonality  relations.  However,  one  can 
obtain  different  "AM" s  depending  upon  the  choice  of  M , 
Also,  the  decision  to  minimize  r,  while  reasonable,  is  never¬ 
theless  arbitrary. 

It  is  evident  that  the  resulting  mass  and  stiffness  matrices 
are  not  unique.  Subsequently,  this  truth  will  he  illustrated  with 
a  simple  numerical  example.  Hence,  it  is  concluded  that  in 
order  to  determine  the  system  matrices  uniquely  ,  some  more 
conditions  in  addition  to  the  (necessary  hut  not  sufficient)  or¬ 
thogonality  conditions  must  be  satislied.  These  additional 
conditions  can  be  readily  obtained  from  the  equations  of  mo¬ 
tion.  The  frec-response  data  can  be  used  to  determine  the 
eigenvalues  and  eigenvectors  of  the  system.  The  estimated 
modal  data  then  can  be  used  in  conjunction  with  the  forced - 
response  data  to  uniquely  identify  the  system  matrices. 


Identification  of  Fjigenvalucs  and  Eigenvectors: 
Rajaram  and  Junkins*  Approach 

A  system  described  by  Eq.  (33)  will  he  considered.  The 
modal  coordinate  transformation  is  introduced  as 


xU )  =  Etiil) 


where  vj ( r )  is  the  normal  or  modal  coordinates  oT  the  system. 
Introducing  Fq.  (39)  into  the  equations  of  motion,  Fq  (33). 


Ml  r)(t )  r  KE\ |(r)  -  0 


Multiplying  Fq  (40)  by  I  ' ,  wc  get 


E' MEiid)  4  E'KEtiU)  0 


Due  to  the  orthogonality  properties  of  the  eigenvectors.  Fq 
(41)  represents  a  set  of  "n"  uncoupled  second-order  cqua 
lions.  It  the  eigenveetors  arc  normalized  with  respect  to  the 
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mass  matrix  as  per  Eqs.  (34),  Eq.  (41)  becomes 

4(0  +  |ur)*(0  =0.  (uj- ]  =diag(u') . w;’)  (42) 

When  the  eigenvectors  are  not  normalized,  Eq.  (41)  can  be 
written  as 

M„,iHO  +  Amij(D=0  (43) 


Thus,  the  parameter  matrices  can  be  estimated  uniquely.  We 
need  to  estimate  only  “n”  parameters,  viz.,  the  diagonal 
elements  of  Mm.  The  elements  of  K,„  can  be  derived  from 
those  of  Mm  through  Eq.  (44).  Also,  the  amount  of  forced- 
response  data  required  to  estimate  the  modal  matrices  is  not 
large.  We  need  only  “2 n”  measurements  ("n"  accelerations 
and  "n"  displacements),  in  addition  to  the  measurement  of 
forces. 


M,„  and  Am  are  diagonal  matrices;  the  generalized  "modal 
mass”  and  generalized  "modal  stiffness  matrix,”  respectively. 
Also,  Am  is  related  to  M,„  by  the  following  relationship; 

M„  =  l  u>- 1  K,„  (44) 

The  solution  of  Eq.  (42)  can  be  written  as 

ti,  (O  =  r,cosu),/  +  j,sinu),f  i=l.2 . n  (45) 

f,  and  s,  are  constants  depending  upon  ij,(0)  and  n,(0). 
Substituting  Eq.  (45)  into  the  transformation  Eq.  (39),  we 
obtain 


x(/)=  (,4,cosu),/ +  S.sinu),/)  (46) 

i  / 


where  it  is  evident 


A  ,=c,e,  and  B,  =  s,e, 


Identifying  either  A,  or  B,  is  equivalent  to  identifying  a  scaled 
version  of  the  /th  normalized  eigenvector.  A  Gauss-Newton 
least -squares  differential  correction  method  or  a  direct 
method  based  on  the  Fourier  transform^  of  x(l)  can  be  used 
to  obtain  the  modal  parameters  {ui,,A,.Bl ' . 

We  now  turn  attention  to  estimating  the  properly  scaled 
mass  and  stiffness  matrices.  Equation  (43),  in  the  presence  of 
forces,  becomes 


Identification  of  Damped  Systems 

We  now  turn  our  attention  to  the  necessary  modifications  of 
the  preceding  approach  for  including  viscous  (or  equivalent 
viscous)  damping.  The  equations  of  motion  for  a  damped 
system  are  given  by  Eq.  (I). 

The  eigenvalues  and  eigenvectors  of  a  damped  system  are 
complex  quantities.  In  order  to  apply  classical  modal  analysis 
techniques,  it  is  a  usual  practice  to  assume  that  the  damping  is 
either  small  or  of  a  proportional  type.  Since  the  measured 
modes  are  complex,  methods  have  been  proposed  to  extract 
the  normal  modes  from  the  complex  modes.  However,  it  is 
possible  to  rigorously  apply  a  generalized  modal  analysis 
technique  by  transforming  Eq.  (1)  from  configuration  to  state 
space  and  estimate  the  system  matrices,  analogous  to  the 
previous  section. 

Introducing  the  state  vector 

g(l)  =  |x'  ( z )jt ( / )  ]  ' 

Equation  (1)  can  be  written  as 

M’gU)  +  *’*(/)  =/•</)  (51) 


where 


-  K  0 
0  M 


(52a) 


A* 


0 

A 


A 

C 


(52b) 


M„lV)  +  Am*(0  =  E’fU)  (47) 

where  the  (nxl)  force  vector  /(t)  may  contain  zero  entries, 
i.e.,  all  of  the  degrees  of  freedom  need  not  be  forced,  M„.  and 
A'„,  are  easily  determined  from  the  scalar  components  of  Eq. 
(47),  using  the  fact  that  Mm  (/,/)  =  wfA'm  (/,/>.  It  should  be 
noted  that,  henceforth,  the  notation  £  will  be  used  to  repre¬ 
sent  the  measured  eigenvectors,  normalized  with  respect  to  the 
a  priori  mass  matrix.  Since  E  is  measured,  transformation 
equation  (39)  can  be  used  to  transform  measurements  in 
physical  space  to  modal  space,  i.e.. 


11 

'X(I) 

<48a) 

y>U)-E 

'xU) 

(48b) 

Introducing  Eqs.  (48)  into  Eq.  (47),  for  a  known  force  vector, 
it  is  obviously  possible  to  determine  the  diagonal  elements  of 
the  modal  stiffness  matrix  A'„,  and,  using  Lq.  (44),  ,Vf,„  can  be 
computed.  The  properly  scaled,  configuration  space-mass 
matrix  then  can  be  obtained  from 

M  ~  E  1  M,„ E  '  (49) 

Similarly  the  stiffness  matrix  is  given  as 

A =£  'A „,£  '  (50) 


r  o  ) 

r'\  /  J 

The  eigenvalues  and  eigenvectors  of  a  system  described  by  Eq. 
(51)  occur  in  complex  conjugate  pairs,  i.e.,  if  A,  is  an  eigen¬ 
value,  A,  is  also  an  eigenvalue.  Similarly  e,  and  e,  are  the 
eigenvectors  of  the  system.  The  orthogonality  relations  arc 

EJ  M'E  =  /  (53a) 

E’K'E-  —  A  (53b) 

where  A  is  a  diagonal  matrix  of  the  eigenvalues.  Note  that  the 


modal  matrix  is  of  order  (2 n  x  2 n).  Also,  the  eigenvectors  arc 
normalized  with  respect  to  M'  to  satisfy 

i  =  1,2 . 2n  (54) 

When  the  eigenvectors  arc  not  normalized,  of  course,  Eqs. 
(53)  become 

E’M'EM:,  (55a) 

E'K'E  A’  (55b) 


For  high-dimensioned  systems,  of  course,  the  inverses  shown  where  M’„,  and  K’n,  arc  diagonal  bin  complex  matrices.  The 

are  replaced  by  appropriate  matrix  reduction  algorithms.  same  notations  as  in  the  previous  section  arc  used.  The 
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eigenvectors  have  the  form 

-U. } 

The  free  response  of  the  system  can  be  written  as 

x(f)  =  Yj  (fl/e*'1'  +  a,ex'')  (57) 

i  - 1 

where  A,(  =  -a,+ju\)  and  the  /th  eigenvalue,  and  a,  is  the 
damping  factor  and  u\  the  damped  frequency  of  oscillation. 
Equation  (57)  also  can  be  written  as 

n 

x (t)=2  ^  e*"i'(C,cosu;,f  —  5,sinu>,r)  (58) 

i  -  i 


Although  Berman’s  corrections  to  the  diagonal  elements  are  in 
the  right  direction,  the  off-diagonal  elements’  corrections  are 
of  comparable  size  and  are  no  longer  zero.  A  significantly  dij- 
ferent  final  mass  matrix  estimate  would  have  been  obtained, 
of  course,  if  M ^  were  chosen  differently.  The  estimated  stiff¬ 
ness  matrix  for  Berman’s  approach  is  found  to  be 

68.36  -32.52 

A  = 

-  32.52  71.68 

It  can  be  noted  that  the  diagonal  terms  are  corrected  fairly  well 
while  the  corrections  to  off-diagonal  terms  are  relatively  small 
in  this  case. 

Using  the  method  developed  above  (with  the  eigenvalues 
and  eigenvectors  calculated  using  a  finite  Fourier  transform), 
the  mass  matrix  is  determined  from  Eq.  (49)  to  be 


C,  and  S,  are  the  real  and  imaginary  components  of  a,,  respec¬ 
tively.  A  Gauss-Newton  least-squares  differential  correction 
method*  can  be  used  to  identify  S,,  a,,  and  u.-, 

(I  =  1,2 . n).  Fast  Fourier  transform  of  x(l)  is  quite  useful  in 

this  case.  The  frequencies  can  be  estimated  from  the  power 
spectral  density  (psd)  plot  and  used  as  a  priori  values  in  the 
Gauss-Newton  algorithm.  In  this  way,  the  convergence  do¬ 
main  of  the  algorithm  can  be  enhanced  considerably  .  Using 
the  orthogonality  relations  [Eqs.  (56)]  and  the  trans¬ 
formation 

gU)  =  [*,T|,(r) +e,r)  (0)  =£ij(r)  (59) 

i  - 1 

Eq.  (52)  reduces  to 

=  A'„t)(r)  +  E'f  U)  (60) 

where 

V(0=  [vill)iil{t)...r,„U)fiAO\’  (61) 


[  100.00  -0.8415E-  04 

-  0.84I5E-  04  200.0 

and  the  estimated  stiffness  matrix,  from  Eq.  (51),  is  calculated 
as 

71.995  -  35.997  ] 

A  = 

-35.997  71.9996  j 

The  small  residual  errors  art  the  consequence  of  truncation  of 
the  Fourier  transform  of  -r(/)  to  obtain  eigensalues  or 
eigenvectors.  If  the  Gauss-Newton  iteration  is  used  instead, 
the  .1/  and  A  matrix  are  recovered  exactly  (to  eight  digits).  It  is 
evident  from  this  simple  example  that  the  proposed  scheme 
correctly  identifies  the  system  matrices  to  within  truncation  er¬ 
rors  in  the  finite  Fourier  transform.  In  essence,  the  scaling  im¬ 
plicit  within  Berman's  correction  norm  minimization  is  re¬ 
placed  by  the  requirement  that  the  estimated  M  and  A  be  con¬ 
sistent  with  a  measured  forced  response. 


is  a  complex  modal  coordinate  sector.  Equation  (60)  can  be 
used  to  identify  the  M’„  and  A‘  matrices  from  forced 
response,  analogous  to  the  previous  section.  After  determin¬ 
ing  the  modal  matrices.  M’  and  A’  can  be  obtained  from  Eqs. 
(56).  The  method  is  similar  to  the  one  for  the  undamped 
system,  except  that  the  quantities  involved  are  complex. 

Numerical  Examples  Using  Free 
and  Forced  Response 

Example  1 

Consider  a  two  mass-spring  system.  The  mass  and  siiflncsx 
are  given  as 


Example  2 

A  two  mass-spring  damper  system  is  considered.  The 
various  matrices  are 

I  I  0  1  15-4]  0.4  -  ().: 

M--  j  .  A  I,  C- 

0  I  -44]  -  0.2  0.2 

The  eigenvalues  are 

X,.X,=  -0.222593  ±y2. 578255 

A...  X ,  -  0.027406  ± yO. 545796 


100  ol  72-36 

M  I  A 

0  2(H)  -  36  72 


Choosing  the  a  priori  mass  matrix  A/.,  and  stiffness  matrix  A  , 
for  Berman's  method  '  1  as 


[90  0|  65  -  32  I 

M.=|  |  A,, 

[  0  220  j  -  32  79  | 

The  true  system  eigenvalues  and  eigenvectors  arc  used  as 
measurements  Alter  carrying  out  algebra  of  Hetman’s 
method,  the  estimated  mass  matrix  is  found  to  be 


[  96.67  6.67 

M- 

6.67  206.67 


Using  the  Gauss-Newton  method,  the  eigenvalues  and 
eigenvectors  arc  estimated  from  free  response.  The  free- 
response  data  are  properly  scaled  by  applying  an  impulsive 
force  on  the  second  mass.  The  system  matrices  are  obtained 
using  the  excitation/,  -  0.1  sign  (sin  0.2r).  The  mass,  stil  lness, 
and  damping  matrices  are  identified  exactly  (eight  digits). 
Thus  the  present  method  generalizes  fully  to  include  arbitrary 
viscous  damping 

Conclusions 

Three  novel  schemes  arc  proposed  to  identify  the 
parameters  ol  vibrating  structures.  Numerical  results  on  a 
vaiietv  of  transparent  examples  support  the  validity  of  all 
three  methods  The  physical  properties  of  mass,  stillness,  and 
damping  matrices  are  identified.  All  three  proposed  methods 
are  applicable  to  damped  structures.  No  assumptions  rcg.it 
ding  the  nature  ol  damping  are  made,  other  than  it  is  of  the 
viscous  type  Systems  with  closely  spaced  frequencies  present 
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no  apparent  computational  difficulty,  in  fact,  example  struc¬ 
tures  with  repeated  frequencies  and  rigid-body  modes  are 
identified  reliably  without  difficulty.  It  is  also  shown  that 
multiple  excitation  vectors  should  be  chosen  to  form  an  in¬ 
dependent  set.  The  methods  have  been  illustrated,  however, 
only  for  low -dimensioned  examples;  significant  future  effort 
should  consider  high-degree-of-freedom  systems  to  evaluate 
the  robustness  and  relative  merits  of  these  approaches.  It  is 
also  important  to  reduce  the  dimensionality  by  coupling  the 
estimation  algorithms  to  the  structural  modeling  process. 
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An  eigenspact  optimization  approach  is  proposed  and  demonstrated  for  the  design  of  feedback  controllers  for 
the  maneuvers  vibration  arrests  of  flexible  structures.  The  algorithm  developed  is  shown  to  be  equally  useful  in 
sequential  or  simultaneous  design  iterations  that  modif»  the  structural  parameters,  sensor/actuator  locations, 
and  control  feedback  gains.  The  approach  is  demonstrated  using  a  differential  equation  model  for  the 
"Draper  RPL  configuration."  This  model  corresponds  to  the  hardware  used  for  experimental  verification  of 
large  flexible  spacecraft  maneuver  controls.  A  number  of  sensor/actuator  configurations  are  studied  vis-a-vis 
the  degree  of  controllability.  Linear  output  feedback  gains  are  determined  using  a  novel  optimization  strategy. 
The  feasibility  of  the  approach  is  established,  but  more  research  and  numerical  studies  are  required  to  extend 
these  ideas  to  truly  high-dimensioned  systems. 


Parameterization  of  the  Controlled  System’s 
Eigenvalues  and  Eigenvectors 

CONSIDER  a  linear  structure  (modeled  by  a  finite  element 
or  similar  discretization  scheme)  in  which  the  configura¬ 
tion  vector  x  is  governed  by  the  system  of  differential  equations 


For  linear  output  feedback,  we  seek  the  constant  gain  matrices 
G,.  G2.  and  G,  so  that 


«  =  -  [ G,y,  +  G,y,  +  G,y,  1 
=  -GiH/X-  G,H2x  -  G,H,x 


(3) 


\tx  +  Or  +  Kx  =  Bu 


(1) 


where 


Substitution  of  Eq.  (3)  into  Eq.  (1)  gives  the  closed-loop 
system 


M  =  n  v  n  symmetric  positive  definite  mass  matrix 
C  -  n  x  n  symmetric  positive  semidefimte  structural  damp¬ 
ing  matrix 

K-nxn  symmetric  positive  semidefinite  stiffness  matrix 
B  =  n  x  m  control  influence  matrix 
x  =  n  x  I  configuration  vector 
u  =  m  x  I  control  vector  (  )  =d/df(  ) 

Considering  the  case  of  linear  output  feedback  control,  let  the 
locJ  position,  velocity,  and  acceleration  measurements  be 
denoted  by 


Mx  +  Cx  =  Kx  =  0 

where  the  closed-loop  system’s  matrices  are 

M=M+BG,H„  C=C+BG:H \,  K  =  K*BG,H, 
Introduce  the  notations 
M  =  Mia).  C=C(o).  K  =  K(a) 

B  =  B{c).  H,  =  H,{b).  G,=C,(g) 


(4) 


(5) 


(6) 


y,-H,x.  y.  =  H,x,  y,  =  H,x 


(2) 


which  represent  the  linear  relationship  of  the  locally  measured 
position  velocity  y  ..  and  acceleration  yt,  where 


y,  -  m,  x  vector 


H ,  =  m ,  x  n  matrix 


y.-m-x  I  vector  H.  -  m,  x  n  matrix 
y,  --  m.  -  1  vector  H ,  =  m,  x  n  matrix 


where  a  is  a  vector  of  the  structural  and  geometric  model 
parameters  (defining  mass,  stiffness,  damping,  configuration 
lengths,  cross-sectional  areas,  etc.),  b  a  vector  of  the  sensot- 
type  and  location  parameters,  c  a  vector  of  the  actuator-ty  pe 
and  location  parameters,  and  g  a  vector  of  the  control  gains 
Defining  the  .V*  I  global  structural  and  control  parameter 
vector  as 


p’  =  \aTbTcrgT\ 

it  is  apparent  from  Eqs.  (5-7)  that 

M=M(p),  C=C(p),  K  =  K(p) 


(7) 


(8) 
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The  details  of  the  parameterizations  of  Eqs  (8)  are  dependent 
upon  the  particular  modeling  approach.  Often,  these  are  sim 
pic  algebraic  expressions  (the  elements  of  p  appear  explicitlv  in 
M,  C,  K) 

Considering  the  first-order  state-space  differential  cqua 
tion,  which  is  the  equivalent  to  the  second-order  closed  loop 
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where  the  conventional  normalizations  for  the  eigenvectors 
are  adopted  as 

tjA0,=&u.  *fB0,  =  6,,  A,  (13) 


dominate  a  measure  of  optimality  or  sensitivity.  This  allows 
judicious  and  lower-dimensional  suboptimizations  to  be  made 
and  is  a  key  to  attacking  truly  high-order  systems. 

These  two  qualitative  ideas  serve  as  the  main  motivation  of 
our  approach.  Of  course,  an  initial  design  point  (values  for  the 
elements  of  p)  is  required.  Typically,  the  initial  design  point 
will  be  the  output  of  some  (arbitrary)  structural  design  process 
and,  probably,  a  constant-gain  optimal  regulator  design  for 
the  given  structure.  Thus,  the  present  family  of  algorithms  are 
designed  to  begin  with  the  typical  output  of  a  conventional  se¬ 
quential  structure — controller  design  process.  However,  for 
moderate-dimensional  applications,  we  have  been  able  to  in¬ 
itiate  control  gain  optimizations  with  a  free-vibration  (zero 
control  gains)  case  and  still  achieve  reliable  convergence. 

Central  to  the  application  of  the  optimization  algorithms 
developed  below  lies  the  necessity  to  compute  efficiently  the 
partial  derivatives  of  the  generalized  eigenvalues  and  eigenvec¬ 
tors  of  Eqs.  (12-14).  Attractive  algebraic  equations  have  been 
derived  that  explicitly  determine  these  derivatives  as  "side 
calculations"  at  much  less  computational  expense  than  the 
solution  of  the  eigenvalue  problem  itself.  The  development  of 
these  equations  is  briefly  summarized  below  for  the  first  and 
second  partial  derivatives  oi  the  eigenvalues. 


Since  A  =A(p)  and  B  =  B(p),  it  seems  natural  to  consider  the 

eigenvalues  {  \ . ^2n  l  and  eigenvectors  1 0,,^, . 02n,'h„  I 

to  be  functions  of  the  parameter  vector  p,  viz., 

X,  =  X,(p).  0,=0,(p),  iL  =  ^,(P)  (14) 

Except  for  occasional  singular  events  (e.g.,  multiple  eigen¬ 
values  or  near-multiple  eigenvalues),  the  nonlinear  functional 
dependence  of  Eqs.  (14)  can  be  assumed  to  be  continuous. 

Qualitative  Approach  to  Eigenspace  Optimization 

Most  structural  and  control  optimality  criteria  can  be  stated 
explicitly  in  terms  of  X,  and  or  directly  as  functions  of  p.  It 
is  obvious  that  an  algorithm  that  can  effectively  optimize  p 
(over  some  admissible  set  to  minimize  some  optimality  criteria 
and  satisfy  the  constraints  stated  as  functions  of  p  and  X,)  pro¬ 
vides  a  direct  method  for  controller/structure  optimization 
problems.  Unfortunately,  there  are  several  formidable  dif¬ 
ficulties,  the  two  most  prominent  being: 

1) The  Nx\  p  vector  is  of  high  dimension,  even  for  struc¬ 
tural  and  control  system  models  of  moderate  complexity.  The 
dimension  N  of  p  can  be  several  hundred. 

2)  The  functional  relationship  implied  by  Eq.  (14)  is  the 
solution  fo  the  large  eigenvalue  problem  of  Eqs.  (12)  and  (13). 
It  is  typically  a  highly  nonlinear  function  of  p  and  is  occa¬ 
sionally  characterized  by  singular  local  behavior  (bifurcation 
points  at  repeated  roots,  for  example). 

Practical  optimization  algorithms  that  can  deal  with  these 
two  sources  of  difficulty  in  a  rigorous  and  globally  convergent 
way  do  not  exist.  However,  we  have  developed  a  strategy  for 
carrying  out  optimizations  and  suboptimizations,  in  spite  of 
these  two  sources  of  difficulty.  There  are  two  heuristic  ideas 
underlying  our  approach: 

1 )  Regions  of  extremely  high  sensitivity  (to  the  p  vector)  are 
generally  undesirable.  Therefore,  if  the  performance  index  or 
constraints  include  a  measure  of  eigenvalue  solution  sensitiv¬ 
ity,  successful  suboptimizations  are  obviously  less  likely  to  en¬ 
counter  the  singular  events.2 

2)  Exact  eigenvalue  placement  (or  “pole  placement”),  for  a 
high-order  system,  is  not  a  reasonable  design  approach. 
Rather  than  attempting  to  prescribe  an  exact  point  location 
for  every  eigenvalue,  it  is  more  reasonable  (and  leads  to  better 
algorithms)  to  move  all  of  the  eigenvalues  into  an  acceptable 
region  of  the  complex  plane.  From  the  viewpoint  of  pole 
placement,  this  allows  attention  to  be  locally  concentrated 
upon  just  the  "problem  children"  eigenvalues  that  are  the  far¬ 
thest  outside  the  acceptable  region  or  those  that  locally 


Partial  Derivatives  of  the  Closed-l.oop  Eigenvalues 
with  Respect  to  Structural  Model 
and  Control  System  Parameters 

Differentiating  Eqs.  (12)  with  respect  to  a  typical  element  p, 
of  p,  upon  premultiplying  the  resulting  two  equations  by  4.J 
and  0j  and  making  use  of  Eqs.  (13),  gives 


d<t>  ,  T  SB  dA  1 
K--K,H,A—^,  [  -K 


r  ,  dt,  , r as'  a,4'i 

X, -\  )4>r/t  — —  =  0  — - X— — U, 

dp,  ■  L  dp,  dp,  J 


X,^k 

dp,  J 

(16) 

,2zr;  y  =  1,2,. 

..,2  n\ 

and  f=  1,2,...,/V.  For  i  =  j,  both  Eqs.  (15)  and  (16)  reduce  to 
the  following  well-known  result1  for  the  gradients  of  the  2 n 


eigenvalues 


dp,  L  dp,  dp,  J 


(17) 


Thus,  having  solved  for  eigenvalues  and  eigenvectors,  a 
moderate  side  calculation  produces  the  first-order  eigenvalue 
sensitivities.  Differentiating  Eq.  (17),  with  respect  to  p^.  we 
obtain  the  following  equation  for  the  second  partial 
derivatives: 


d-'\  d'B  ax,  8.4  8-4  1 

- - 1_  =  iaH - - - X, -  o, 

dp,dp„,  '  l  dp,dpm  dp„,  dp,  dp,dp„,  J 


i^[^--xil]*+W— -a-1^ 

dpm  L  dp,  dp,  J  L  dp,  dp,  J  dp„, 


(18) 


Since  Eq.  (18)  involves  d0,/dp„  and  di,'dp„,  we  must  either 
evaluate  or  eliminate  them,  Wc  choose  the  latter  method 
Following  Plaut,1  we  project  8^,/8/im  and  d0,'dpm  onto  the 
eigenvectors  themselves  as 


d0, 

dp, 


2n 

Y  “».*«• 

*  -  i 


H, 

dp, 


Y  btAi 


(19) 
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where  a,,  and  bk,  are  scalar  constants.  For  ixj,  substituting 
Eqs.  (19)  into  Eqs.  (15)  and  (16)  yields 


For  the  case  of  distinct  eigenvalues,  Eqs.  (20)  and  (21)  provide 
the  a,,  and  bp  values,  except  for  a„  and  b„  (which  remain 
undetermined).  The  normalization  of  Eq.  (13),  upon  setting 
y  =  /and  differentiating  with  respect  top,  and  substituting  Eq. 
(19),  yields 


+  b„  =  -  4J- 


wmll  '  ''It  T  I  a  ~ I 

dp, 

Substitution  of  Eq.  (17)  into  Eq.  (18)  gives 
3-X,  d'B  d- A  1 


3-X,  ,  f  d'B  ,  d- A  1 

dp,dp„  '  L  dp,dp„,  ap(dpm  J  ' 

,  dA  ( ai,  3  t  r  a<t>, 

-  4!G,„,<t>,4  o,  +  j  —  (23) 

dp,  C  dpm  )  dpm 


with  G,(  ^dB/dpk  -\,(dA/dpt).  Substitution  of  Eqs.  (19) 
into  Eq.  (23)  gives 

3’X,  .  ,f  d  B  d:  A  I  .  _  „  dA 


a-'K  ,  r  a  b  a-' a  l 

=  i!  -z — 7 - x,— —  ♦.  - 

L  dp,ap„,  ap,ap„  J 


+  («„  +  b„)4?G„o,  +  ldt,4[G„<t>,  +  at,4,rG,t4>  J 


Finally,  eliminating  a„  +  using  Eq.  (22)  and  a,,  and  b;,  for 
jxi  using  Eqs.  (20)  and  (21 )  gives  the  final  expressions  for  the 
second  partial  derivatives, 

a- x,  ,  r  d-B  a- a  1 

- —  =  j,n - x, - —  o, 

dp,ap„,  L  dp,dp„,  dp,dpm  J 


,dA  r  dA 

-  4!  g,„,0,4! — — o,  —  4!  -<u,VG,,a> 

ap,  a  pm 


t  [* 

i  f.i  xt  1 


G„„<t>k  4l G„<t>,  +  VG,„,<i,V G„<t>k 


1  <  X,  -  x* )  J 

which,  of  course,  are  singular  for  repeated  eigenvalues. 


Closed-Loop  Response 

The  response  of  the  controlled  configuation  to  initial  distur¬ 
bances  is  determined  using  a  modal  approach.  The  modal 
coordinates  ij  are  mapped  into  the  state  space  by 


Z=l©]>),  Z  =  1 0 1 1» 


where  (0j  is  the  right  modal  matrix.  Substituting  Eqs.  (26) 
into  Eq.  (9),  premultiplying  by  (i£)r  (the  left  modal  matrix) 
and  utilizing  Eqs.  ( 1 3),  the  equations  of  motion  are  uncoupled 
and  given  by 


V,  =\v,. 

The  solution  to  Eq.  (27)  is 
y,=eK<"  '«'ij,(r„), 


/  =  1,2 . 2n 


i-  l,2,...2n 


The  initial  conditions  are  mapped  into  modal  space  by 
premultiplying  Eqs.  (26)  by  |0)  1  =  |^]  ’A  and  utilizing  Eq. 
(13)  to  give 


y(io)  =  WTAz(i0) 


The  system  response  in  normal  coordinates  is  then  given  by 
Eqs.  (28)  and  (29).  However,  the  eigenvalues  and  eigenvectors 
are  in  the  general  complex  consisting  of  n  eigenvalues  and  n 
eigenvectors  plus  their  associated  n  complex  conjugates.  We 
seek  a  solution  that  will  eliminate  the  complex  conjugates  and 
thus  allow  us  to  truncate  the  number  of  modes  used  in  the 
solution. 

The  modal  matrices  are  partitioned  as  follows: 


' 

l,  • 

<t>l 

[V  = 

.  1  0  1  = 

t h 

1}  - 

0, 

*2  . 

where  4,.  4  0(,  and  <t>2  are  nxk  matrices  normalized  by 

Eqs.  (13).  It  follows  that  4i,  <f >,.  and  0.,  are  also  nxk 
matrices  and  are  thus  the  complex  conjugates  of  4,.  4 0.’, 
respectively.  The  vector  of  normal  coordinates  y  can  also  be 
partitioned  as 

'*{  t } 


where  f  is  a  A:  x  1  vector.  It  can  be  shown  that  f  is  also  k  x  1 
and  is  the  complex  conjugate  of  f.  Substituting  Eqs.  (30)  and 
(31)  into  Eq.  (26),  the  response  in  configuration  coordinates  is 


f  *  1  f  0/r+0,f 

Z  =  1  !• =  1  ■ 
L.  X  J  l  0,f  +  0,f 


^ f  I Be<j>i)  |/fef|  -  [/m0,]  |/mf|  j 

1  [Red;]  I /fefl  -  [lm<t>;]  |/mf|  J 


From  Eq.  (32),  it  is  evident  that  the  complex  conjugates  of  the 
n  eigenvalues  and  eigenvectors  are  not  needed  and  that  the  k 
modes  can  be  used  to  determine  the  response.  The  measured 
beam  deflections  and  deflection  rates  follow  from  Eq.  (2)  and 
the  controlled  response  from  Eq.  (3). 


A  Model  of  Draper/RPL  Configuration 

Referring  to  Figs.  I  and  2,  we  consider  the  planar  rota¬ 
tional/vibrational  dynamics  of  the  demonstration  experimen¬ 
tal  model  sponsored  by  the  U.S.  Air  Force  Rocket  Propulsion 
Laboratory  for  testing  the  control  laws  for  maneuvering  flexi¬ 
ble  spacecraft;  this  experimental  work  is  presently  being  con¬ 
ducted  at  the  C.S.  Draper  Laboratory. :  The  central  hub  is 
supported  by  an  air  bearing  table  with  four  appendages  can¬ 
tilevered  from  the  hub.  Table  1  summarizes  the  Draper/RPL 
configuration  parameters.  Following  Turner  and  Chun,:  we 
form  a  discretized  model  for  the  system  by  assuming  that  the 
elastic  deformations  of  each  of  the  arms  (relative  to  a  bodv- 
fixed  undeformed  state)  can  be  represented  as  linear  combina¬ 
tion  of  the  comparison  functions. 


0A<z)  =  /-cos(^p)  +-L< -/)*•' 


so  that  the  transverse  body-fixed  deformation  of  the yth  arm  is 
modeled  as 


y,(f.Z)  =  E<M')*»(I).  j~  2.2, 3, 4,  0<;</  (34) 


yy 
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Radial  elongation  of  the  arms  is  neglected,  as  are  out-of-plane 
deformations.  For  the  numerical  results  presented  below,  we 
took  /V=  10.  Using  Eq.  (34)  to  evaluate  the  potential  and 
kinetic  energy’  leads  to  equations  of  the  form  of  Eq.  (1),  with 
the  configuration  vector 

4=  [6 \Q  nQ  2i""  :Q\i \QnQ2i""Qs2  \'“'Q  nQ  n'"Q  S4 1 1 

We  restrict  attention  to  the  class  of  antisymmetric  deforma¬ 
tions  whereby  y,U,z)"y2(t.z)  and  y3{i.z)  *>,(f.z);  thus, 
Q,iU)mQ,2(i )  and  ql}(t)*ql4(t).  For  this  class  of  antisym¬ 
metric  motions,  the  configuration  vector  can  be  collapsed  to 

4=  [0\q,,q2i-Qsi\QiiQ23- -Qs)\J  <35) 

For  N=  10,  the  order  of  the  system  of  Eq.  (1)  is  thus  21.  The 
explicit  expressions  for  the  elements  of  the  M  and  K  matrices 
are  developed  in  Ref.  3.  We  take  C-(iK,  where  n  is  a  con¬ 
stant.  For  the  numerical  examples  presented,  p  =  1  •£-  5;  the 
controlled  response  of  the  first  seven  modes  is  insensitive  to  an 
order  of  magnitude  variation  of  p.  In  Fig.  3,  we  show  the  first 
nine  normal  modes. 

The  eigenvalues  and  eigenvectors  for  the  first  seven  modes 
have  converged  in  the  numerical  sense  that  increasing  N  does 
not  change  the  first  four  or  five  digits,  whereas  modes  eight 
and  nine  have  converged  to  about  three  digits.  (However,  the 
last  ten  higher-frequency  computed  eigenvalues  and  eigenvec¬ 


tors,  as  expected,  are  not  as  accurately  calculated.)  Thus,  we 
restrict  our  optimization  discussions  to  the  first  nine  modes. 
After  the  first  (rigid-body)  mode,  the  remaining  modes  occur 
in  near  pairs.  The  antisymmetic,  “opposition”  modes  are  sim¬ 
ple  cantilever  beam  modes  characterized  by  the  adjacent 
beams  moving  in  opposition  (the  constraint  torques  between 
the  hub  and  the  beam  clamp  cancel  in  equal  and  opposite 
pairs,  resulting  in  a  zero  rotation  of  the  hub).  The  antisym¬ 
metric  “unison”  modes  are  perturbed  cantilever  modes  (with 
varied  frequencies)  characterized  by  all  four  beams  moving  in 
unison;  the  hub  has  nonzero  rotation  for  these  modes.  As  is 
evident,  the  corresponding  pair  of  unison  and  opposition  fre¬ 
quencies  are  closely  spaced;  this  spacing  decreases  with  in¬ 
creasing  mode  number. 


Table  1  Draper/RPL  configuration  parameters 


Hub  radius 

lit 

Rotary  inertia  of  hub.  r 

8  slug-fr 

Mass  density  of  beams,  p 

5.22  slug  ft’ 

Elastic  modulus  of  the  arms,  E 

1.1  x  10  lb  in.2 

Arm  thickness,  ta 

0.125  in. 

Arm  height,  ha 

6  in. 

Arm  length,  L 

4  ft 

Tip  mass 

0.156941  slug 

Rotary  inertia  of  tip  masses 

O.OOIK  slug-ft: 

Fig.  2  A nli-«ym metric  deformation. 


Fig.  3  First  nine  normal  modes. 
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Actuator/Sensor  Configuration 

We  admit  torque  actuators  on  the  central  hub  and  at  some 
station  c  on  each  appendage.  We  also  admit  rotational  posi¬ 
tion  and  velocity  sensors  on  the  hub  and  deflection  and  deflec¬ 
tion  rate  sensors  at  some  stations  on  each  appendage.  For  the 
case  in  which  the  actuators  consist  of  a  torque  u,  applied  to 
the  hub,  a  torque  u,  applied  at  station  z,  on  appendages  I  and 
2,  and  a  torque  us  applied  at  station  z >  on  appendages  3  and  4, 
the  right-hand  side  of  Eq.  (1)  is 


(36) 


\  1 

2 

2  1 

1  1 

1 

Bu  = 

0 

2©  U, 

)  0 

0 

0 

2© ' ( Z  ■ ) . 

l  “r  J 

where 


©  iz)  =  —  [©(c)  I 

dc 


©  U)=  [©,  (z). ..©„,(;)  ] ' 


and,  if  rotational  position  and  velocity  sensors  are  located  on 
the  hub,  while  colocated  deflectiona!  position  and  velocity 

sensors  are  located  at  stations  x, . x „  on  each  appendage 

then  the  sensor  influence  matrices  Ht,  H.  are  both  the  same 
9x21  matrix. 


1  0f  0' 

0  ©r(.v,)  0r 


H,= 


0  ©'<•*,)  O' 

0  O'  ©r(A«) 

0  O'  ©r(*„> 


.  i=1.2 


(37) 


and 


yj  =  1 0y,  (/, X,  )...y,v,x4  )>•,(/, A,  )...>,  (r.Xj, )  ],  y3=  —  (y, ) 


dr 


(38) 


and  the  gains  G,  and  C2  are  both  3x9  matrcies.  For  the 
numerical  examples  below,  we  initially  set  Z/  =z2  =  L/2. 
x,  =xs  =  L/ 4,  x2  =x6  =  L/2,  x3  =x7  =  0.7 L,  x4  =xs  =  0.9 L,  and 
selectively  admit  structural  parameters,  actuator  locations, 
and  sensor  locations  along  with  the  control  gain  vector  in  p. 


Minimum  Modification  Strategy 
for  Structural/Controller  Design  Iterations 

Consider  a  constrained  optimization  problem  wherein  we 
seek  the  optimal  value  of  the  parameter  vector  p  that  ex- 
tremizes  some  performance  measure 


J  =  J\X,(p) . A,„(p),  <t>,(p) . <t>2(p).p] 


(39) 


subject  to  the  satisfaction  of  the  equality  constraints 

or,[A,(p) . A  ;„(p).  ©,<p> . «;„(p).p|  =0 

j=l,2 . N„  (40) 


and  inequality  constraints 

0ji  </3,(A,(p)...,A,„(p),  ©,(p).  .,©>(p),p!  </3,( 

j=1.2 . Nj  (41) 

Thus,  the  performance  measure  and  constraints  are  defined  in 
terms  of  the  eigensolution,  but  we  also  admit  explicit 
dependence  upon  p  to  include,  for  example,  structure  and  con¬ 
trol  system  criteria.  Equations  (39-41)  define  a  nonlinear  pro¬ 
gramming  problem  for  which  a  number  of  algorithms  have 
been  developed  and  applied  during  the  past  two  decades. ** 
One  iteration  stategy  confines  local  attention  to  only  the  lo¬ 
cally  violated  inequality  constraints  and  all  of  the  functions  of 
Eqs.  (39)  and  (40).  Specifically,  one  can  seek  the  smallest  cor¬ 
rection  vector  Ap  that  achieves  specified  increments  of  AJ, 
Aa,,  and  A/3,  for  a  subset  of  the  functions  of  Eqs.  (39-41). 
Linearizing  these  equations  about  p,  results  in 


Ay 


Hrl,]*’ 


(42) 


where  y  =  [•/,<*, ,/3,  ]  T.  Since  Eqs.  (42)  constitute,  typically,  a 
small  number  of  equations  in  a  large  number  of  unknowns,  we 
expect  an  infinity  of  exact  Ap  solutions;  some  criterion  must 
be  introduced  to  select  a  particular  solution.  Motivated  by  the 
desire  to  satisfy,  as  nearly  as  possible,  the  implicit  local  linear¬ 
ity  assumption,  we  seek  a  “small  Ap"  solution.  Minimizing 
the  correction  norm  ApTWAp  (for  IF' a  suitable  weight  matrix) 
subject  to  Eq.  (42)  gives 


Ap  =  W 


not 


’[-tHT1 


(43a) 


Fig.  4  Eigenvalue  placement  locus  of  first  nine  modes:  case  1  (con¬ 
tinuation  parameter  a  varys  from  zero  to  one  in  six  Increments  from 
right  to  left  along  each  eigenvalue  trajectory). 
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where  A  is  a  Lagrange  multiplier  vector  obtained  by  solving 

<4j'” 


Thus,  the  size  of  the  linear  system  we  must  solve  is  equal  to  the 
number  of  functions  [from  Eqs.  (39-41)]  that  we  seek  to 
change  by  the  increment  vector  Ay,  regardless  of  the  number 
of  elements  in  the  A p  vector,  The  partial  derivatives  needed  in 
Eq.  (42)  can  be  evaluated  by  the  chain  rule  partial  differentia¬ 
tion  of  Eqs.  (38-40),  making  use  of  the  eigenvalue  gradients  of 
Eq.  (17).  Of  course,  second-order  optimization  algorithms  can 
also  be  used,5  in  which  case  one  will  need  the  second  partial 
derivatives  of  Eqs.  (25). 

As  a  nonlinear  programming  strategy,  we  first  seek  a 
nearest  feasible  p,  which  satisfies  all  of  the  constraints  [Eqs. 
(40)  and  (41)]  and  then  include  J  increments  in  Ay  to  seek  a 
constrained  optimal  solution.  This  strategy  can  be  formal¬ 
ized4;  it  is  fully  equivalent  to  the  gradient  projection  con¬ 
strained  optimization  algorithms. 

For  example,  we  can  apply  the  above  developments  to  place 
the  Draper/RPL  system's  closed-loop  eigenvalues  in  a  desired 
region  and,  subject  to  this  condition,  minimize  a  robustness 
measure  (e.g.,  the  sensitivity  of  the  eigenvalues  with  respect  to 
variation  of  uncertain  system  parameters). 


Numerical  Examples:  Eigenvalue  Placement 
for  the  Draper/RPL  Configuration 


Case  1 


This  first  example  is  a  modification  of  a  result  presented  in 
Ref.  2.  We  consider  the  problem  of  finding  a  minimum  norm 
gain  vector  gmp  (42  elements  of  G,  and  C2)  that  results  in  the 
eigenvalues  of  the  closed-loop  system  [Eq.  (9)]  satisfying  the 
prescribed  constraints  of 


->(*)  =0  (44) 

where,  in  particular,  we  consider  the  constraint  vector 

><*>=  («>;<*>  !-,<*>  fz<*>-M*)]r  (45) 


The  damping  factors  f,  (g)  and  damped  frequencies  w,  (g)  are 
related  to  the  eigenvalues  X,  (g)  as 


f,=  -tfe[X, (*))/]  [(tfeX, (*))>  +  [/m\  (*)]']  ■ 
ui,  = /m[X,(g)],  i=l,2,... 

The  eigenvalues  are  labeled  according  to  the  ordering 
\Im\,(g)  1  £  ImX^lg)  !<...<  1/znX,  (g)  I 


(46) 


(47) 


The  objective  values  of  the  damping  factors  and  the  first 
natural  frequency  (the  elements  of  y0t,IKmt),  for  the  specific 


•4 


<r 

(0*9) 


16 


tlm  (s) 

Fig.  5  Closed-loop  response  of  Hi):  case  I. 


32 


numerical  example  here,  are  prescribed  as 


7 objective  =  [3.0  0.03  0.03  0.03  0.01  0.01  0.002 
0.002  0.0015  0.0015  ] 


(48) 


In  order  to  enhance  the  convergence,  a  continuation  pro¬ 
cedure  is  used:  Eq.  (44)  is  replaced  by  the  one-parameter  (a) 
family 


<nrobieci.ve-Y  [*(<*>!  =0,  Ozasl 


(49) 


Obviously,  a  =  0  results  in  the  trivial  solution  g(0)  =  0  (cor¬ 
responding  to  free,  uncontrolled  vibration),  whereas  g(  I )  is  the 
desired  solution  [since  Eq.  (49)  becomes  identical  to  Eq. 
(44)] .  Sweeping  a  from  0  to  1  allows  us  to  define  “stepping- 
stone”  problems  that  are  arbitrarily  near  the  converged 
neighboring  solutions;  thus,  the  generalized  Newton 
algorithms  using  Eq.  (43),  with  W'  =  /  and  p  =  g,  can  be  initi¬ 
ated  with  an  arbitrarily  close  starting  estimate  and  very  nearly 
guarantee  satisfaction  of  the  implicit  linearity  assumption.  For 
the  particular  calculations  herein,  we  found  rather  large  a  in¬ 
crements  of  1/6  led  to  reliable  convergence.  Thus,  six  in¬ 
termediate  a  solutions  were  required  to  achieve  final  con¬ 
vergence;  for  each  a  value,  two  or  three  iterations  [Eq.  (43)] 
were  required  to  find  the  g(a)  satisfying  Eq.  (49). 

As  a  specific  example,  we  used  the  initial  g  vector  (displayed 
as  the  elements  of  G,  and  G ,), 


G,=G}  = 


0.001  0  0  0  0  0  0  0  0 
0  00000000 

0  00000000 


The  (1,1)  elements  were  set  slightly  nonzero,  since  exactly  zero 
causes  the  “rigid-body"  eigenvalue  to  be  zero  with  a  resulting 


~  1.2 


X/L .  0.50 

Beam  1  — 
Bean  3  — - 


-1.2  J 


16 


32 


tine  (s) 

Fig.  A  Beam  deflection  response  (midspan  stations):  rase  I. 


1.2 


Ifc  41.  . 

l-iin  ti  An  A  Ax  A  /V'-z- 

;l  i'll  M  1 


IV 


-1.2 


x/l.  i.oo 

Beam  1  - 

Beam  3  — - 


TT 


time  (s) 

Fig.  7  Beam  deflection  response  (tip  stations):  rase 
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Fig.  10  Control  torque  V3U):  cue  I. 


Fig.  0  Control  torque  (.',(():  case  I . 


Fig.  9  Control  torque  l/j(r):  case  I.  Fig.  11  Closed-loop  response  of  HD:  case  2. 


rank-deficient  linear  system  (because  the  damping  ratio  sensitivity  vanishes  if  X,— 0). 
The  final  converged  gain  matrices  were  found  to  be  (only  three  digits  are  shown) 


153 

4.38 

.0  21.5 

.0 

4.38  0.0 

21.5 

0.0] 

0.0 

2.66 

7.89  13.9 

22.5 

2.65  7.89 

13.9 

22.5 

0.0 

2.66 

7.89  13.9 

22.5 

2.65  7.89 

13.9 

22.5 

4.06 

-0.389 

0.0 

-1.73 

0.0 

0.389 

0.0 

-  1.73  0.0 

0.0 

0.094 

-  0.036 

0.188 

-0.100 

0.061 

-0.170 

0.135  -0.211 

0.0 

0.061 

0.170 

0.135 

-0.211 

0.094 

-  0.036 

0.188-0.100 

(50a) 


(50b) 


where  the  12  zeros  indicate  elements  of  G ,  and  G?  not  used  in 
p.  The  locii  of  the  first  nine  closed-loop  eigenvalues  are  plot¬ 
ted  in  Fig.  4  for  0  <oSl,  with  a  =  0  the  point  nearest  the 
imaginary  axis  in  all  cases.  Note  that,  for  eigenvalues  8  and  9, 
the  structural  damping  produced  damping  factors  exceeding 
the  objective  damping  factors  and  thus  these  eigenvalues  ex¬ 
perienced  little  movement.  All  other  higher  calculated  eigen¬ 
values  (Eqs.  (10-12)]  remained  in  the  left  half-plane,  although 
this  does  not  occur  if  the  structural  damping  is  assumed 
negligible. 

The  closed-loop  response  for  the  controlled  configuration 
was  calculated  for  initial  disturbances  of  0(f0)  =5  deg  rigid- 
body  rotation  and  qn(t0)  =0.02,  which  corresponds  to  a  tip 
deflection  of  1.66  in.  for  arms  1  and  2.  The  hub  angle  time 
history  and  beam  deflection  time  histories  for  z/L  =  ‘A  and  1 
are  illustrated  in  Figs.  5-7,  respectively.  The  torque  actuator 
responses  are  shown  in  Figs.  8-10  for  actuators  u,,  u2,  and  u3, 
respectively.  In  Fig.  8,  it  can  be  seen  that  there  are  high- 
frequency  oscillations  superimposed  on  the  first  few  cycles  of 
the  low-frequency  response.  This  is  due  to  the  high-frequency 
vibration  of  modes  8  and  9,  but  these  are  quickly  damped  out 
due  to  the  structural  and  controller  damping.  Also,  note  that 
u,  is  doing  most  of  the  work  compared  to  u2  and  us,  a  desired 


result,  even  though  we  have  not  forced  this  condition  by 
weighting. 

Case  2 

This  case  is  the  same  as  case  I  except  that  we  modified  the 
parameter  vector  to  include  the  sensor  and  actuator  locations, 
in  addition  to  the  42  control  gains  of  case  1,  as 

p  -  actuator  location  vector  c 
=  sensor  location  vector  b 
=  gain  vector  g 

Thus,  the  two  appendage  actuator  stations  ( z,,z2 )  plus  the 

eight  actuator  sensor  ( x . . xt)  stations  brings  the  total 

dimension  of  p  to  52.  Upon  applying  the  minimum  norm  cor¬ 
rection  of  Eq.  (43)  to  satisfy  the  constraint  of  Eq.  (44),  with 
W  =  /,  we  found  the  sensor  and  actuator  stations  were  moving 
undesirably  large  amounts,  so  we  introduced  weights  of  10:  on 
the  actuator  positions  and  the  two  inboard  sensor  stations, 
while  a  weight  of  10s  was  applied  to  the  two  outboard  sensors. 
With  this  modest  artwork  on  the  weight  selection,  the 
minimum  norm  algorithm  was  applied  (with  six  continuation 
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steps  as  in  case  1)  and  reliable  covergenee  ensued  to  place  the  eigenvalues  to  satisfy  Eqs.  (44)  and  (48). 
The  resulting  converged  gain  matrices  were  found  to  be 


G,= 


G-  = 


and  the  converged  (initial  values  in  parenthesis)  appendage  ac¬ 
tuator  stations  were 

Z/  =0.1959  L  (0.5  L),  on  appendages  1  and  2 

z:  =  0.1841  L  (0.5  L),  on  appendages  3  and  4  (53) 

whereas  the  converged  (initial)  sensor  stations  were,  on  appen¬ 
dages  I  and  2, 

x i  =  0.2546  (  0.25  L)  x,  =0.7155  (0.7 L) 

x,  =  0.5414  (0.5  L  )  x4  =0.9317  (0.9  L)  (54a) 

and  on  appendages  3  and  4, 

jr,  =0.2564  (  0.25  L) 
x6  =0.5443  (  0.5  L) 
x7  =  0.7167  (0.7  L) 

x,  =  0.9321  (0.9  L)  (54b) 
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0.165 

0.0 

-2.43 

0.0 

-0.305 

0.0 

-2.80 

0.0  j 

0.0 

- 

1.09 
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-2.80 

4.54 

-  1.13 

1.43 
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3.87 

0.0 

- 

1.10 

1.85 

-2.60 

4.28 

-  1.21 

1.72 

-3.04 

4.46 

(52b) 

The  continuation/root  locus  of  the  eigenvalues  is  essentially 
identical  to  Fig.  4.  However,  the  freedom  to  move  the  sensor 
and  actuator  locations  has  proved  constructive;  by  com¬ 
parison  of  Eqs.  (50)  and  (52).  it  is  obvious  that  the  gains  have 
been  substantially  reduced.  Note  in  Eq.  (53)  that  the  appen¬ 
dage  torquers  hae  been  moved  much  closer  to  the  hub, 
whereas  the  sensors  have  been  modcd  significantly  away  from 
the  hub.  The  net  effect  is  that,  even  though  the  closed-loop 
eigenvalues  have  the  same  position,  smaller  control  torques 
are  required.  This  is  evident  in  comparing  the  controlled 
response  of  case  2  to  the  same  initial  conditions  as  case  I  (Figs. 
11-16),  with  the  corresponding  figures  of  case  1  (Figs.  5-10). 
Note  that  the  peak  torque  is  17  ft-lb  for  case  1,  whereas  it  is 
only  13  ft-lb  for  case  2.  The  controlled  response  is  virtually  un¬ 
changed  (as  might  be  expected,  since  the  two  sets  of  eigen¬ 
values  are  the  same). 

Case  3 

This  case  is  the  same  as  case  2  with  two  new  ingredients; 
I)  two  structural  parameters  the  appendage  length  L  and  the 
tip  mass  mur  are  varied  to  provide  some  structural  design  con¬ 
trol  over  the  free-vibration  frequency  spectrum  of  the  struc¬ 
ture;  and  2)  in  comparison  to  case  2,  the  case  3  controlled 
"rigid-body  mode"  eigenvalue  is  constrained  to  be  the  lower 
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Fig.  18  Beam  deflection  response  (midspan  stations):  cast  3. 


-4. 


0 

Fig.  17 


16 


time  (s) 


32 


Closed-loop  response  of  #(/):  case  3. 


frequency  and  is  much  more  heavily  damped,  to  be  consistent 
with  a  “slewing”  attitude  maneuver/vibration  arrest  control 
law. 

Thus,  we  adopt  the  following  12  objective  constraints  for 
the  closed-loop  eigenvalues  (where  the  nonzero  numbers  in 
parentheses  indicate  initial  natural  frequencies  of  original  un¬ 
controlled  configuration): 


Tobjecme 

'  u,  =  0.3  r/s 
u>2  =  4.5  r/s 
uij  =8.3  r/s 
t,=0. 7 
fz  =  0  03 
<;,  =  0.03 
.  {<  =  0.01 
fj  =  0. 01 
[6  =  0.002 
{,  =  0.002 
{s  =  0.0015 
=  0.0015 


( 0 )  rigid  body  mode  frequency 

(4. 37)  first  flexural  mode  frequency 
(7. 91)  second  flexural  mode  frequency 
(0)  rigid  body  mode  damping  factor 
(0)  first  flexible  mode  damping  factor 
(0)  second  flexible  mode  damping  factor 
(0)  third  flexible  mode  damping  factor 
(0)  fourth  flexible  mode  damping  factor 
(0)  fifth  flexible  mode  damping  factor 
(0)  sixth  flexible  mode  damping  factor 
(0)  seventh  flexible  mode  damping  factor 
(0)  eighth  flexible  mode  damping  factor 


(55) 

These  constraints  were  imposed  in  two  stages.  First,  the  struc¬ 
tural  parameters  (L,MUV )  were  adjusted  to  drive  u2  and  to 
their  objective  values.  In  the  second  stage,  the  52  control  gains 
and  scnsor/actuator  stations  are  modified  (in  6  continuation 
steps  as  in  case  2)  to  drive  y  to  the  above  objective  values. 
Convergence  was  reliably  obtained  and  the  resulting  control 
were  found  to  be 


The  converged  (initial)  appendage  actuator  positions  were 

Zi  = 0.4102  L  ( 0.5  L ),  on  appendages  1  and  2 

Z2  =  0.4102  L  (0.5  L),  on  appendages  3  and  4  (57) 

whereas  the  converged  (initial)  appendage  sensor  positions 
were,  on  appendages  1  and  2, 

x,  =0.2539 L  (0.25  L) 

x2  =  0.5010  L  (0.5  L) 

x,  =0.7608  L  (0.7  L) 

x,  =  0.9653  L  (0.9  L)  (58a) 

and  on  appendages  3  and  4, 

x,  =0.2539  L(0.2S  L ) 

x6  =  0.5010  L  (0.5  L) 

x,  =  0.7608  L  (0.7  L) 

Xg  =  0.9653  L  (0.9  L)  (58b) 
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■  4.84 
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0.0 
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0.037  0.125  0.231 

0.0  ' 

0.341 

(56a) 
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The  converged  (initial)  values  for  L  and  A/lip  were  found  to  be 
L  =  3.67  ft  (4  ft) 


Afnp  =  0.19871  slugs 


(0.156941  slugs) 


The  controlled  response  for  case  3  is  given  in  Figs.  17-22  for 
comparison  with  cases  1  and  2.  Note  the  following  features:  1) 
the  rigid-body  maneuver  and  complete  vibration  arrest  is  ac¬ 
complished  in  approximately  15  s  (about  half  the  time  re¬ 
quired  for  cases  1  and  2);  2)  the  peak  control  torque  is  reduced 
by  about  an  order  of  magnitude;  and  3)  the  location  of  sensors 
is  not  significantly  different  from  case  2,  but  the  actuator  posi¬ 


tions  are  different  (closer  to  midspan).  It  is  also  evident  that 
the  present  maneuver  controls  are  much  more  attractive  than 
cases  I  and  2.  We  have  done  a  variety  of  parameter  variations 
further  establishing  that  the  rigid-body  eigenvalue  placement 
is  the  most  important  feature  of  case  3.  For  any  reasonable 
variation  of  the  sensor/actuator  stations  and/or  structural 
parameters,  we  can  optimize  the  control  gains  to  place  the 
eigenvalues  in  the  same  position  and  achieve  analogous 
results. 


Conclusions 

We  have  developed  and  demonstated  a  minimum  modifica¬ 
tion  strategy  for  structural  and  control  parameter  optimiza¬ 
tion.  Numerical  experience  indicates  occasional  difficulties 
with  the  uncontrolled  modes  being  destabilized,  but  since  all 
of  the  eigenvalues  (up  through  a  conservative  number)  are 
calculated  on  each  iteration,  these  problems  can  be  cir¬ 
cumvented  by  introducing  appropriate  constraints  to  stabilize 
these  modes.  The  method  has  worked  well  on  the  problems 
studied  to  date  and  appears  to  be  an  attractive  approach  to 
follow  in  the  development  of  an  interactive  software  system 
for  structure/controller  design  iterations.  Extension  of  the 
ideas  of  this  paper  to  determine  the  weight  matrices  for  impos¬ 
ing  eigenvalue  placement  upon  optimal  quadratic  regulators  is 
presented  in  Ref.  8.  The  algorithms  presented  herein  have 
been  found  relatively  immune  to  the  high  dimensionality  and 
nonlinearity  of  the  eigenvalue  placement  problem. 
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A  sequential  linear  programming  approach  for 
optimal  placement/constralned  optimization  of 
eigenvalues  and  eigenvectors  of  linear  dynamical 
systems  Is  presented.  As  an  example,  the  total  mass 
of  a  structure  is  minimized  while  the  natural 
frequencies  for  selected  modes  are  gradually  driven 
to  desired  values.  This  highly  nonlinear 
constrained  optimization  problem  is  locally 
linearized  and  linear  assumptions  enforced  by 
specifying  maximum  allowable  local  parameter 
changes.  However,  the  above  constraints  on  the 
magnitude  of  local  parameter  change  restricts  the 
magnitude  of  changes  In  the  system  characteristics, 
and  in  particular,  eigenvalue  constraint  objectives 
which  may  differ  significantly  from  nominal  starting 
values.  The  above  difficulty  is  overcome  by  a  well 
tested  continuation  technique  which  replaces  the 
original,  possibly  rigid,  constraints  by  an 
adjustable  sequence  of  neighboring  constraints.  The 
above  approach  appears  computationally  suitable  for 
redesign  of  high-dimensioned,  complex  dynamical 
systems.  Numerical  examples  are  included  to 
demonstrate  the  practical  merit  of  this  approach. 


NOMENCLATURE 

cost  function 

design  vector 

starting  design  vector 

design  vector  change 

constraint  vector 

constraint  objective  vector 

current  constraint  objective  vector 

continuation  parameter 

maximum  allowable  local  parameter 

change  vector 

transformed  design  vector  change 

generalized  displacement  vector 

global  mass  matrix 

global  stiffness  matrix 

mass  matrix  for  J-th  element 

stiffness  matrix  for  J-th  element 

beam  element  length 

mass/length  distribution 

bending  stiffness  distribution 

Young's  modulus 

mass  density 

element  shape  function 

nondlmenslonallzed  element  coordinates 

lower  bound  vector  on  nodal  thickness 

upper  bound  vector  on  nodal  thickness 

eigenvalue  vector 

objective  eigenvalue  vector 

natural  frequency 

eigenvector 
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INTRODUCTION 

The  practical  need  for  optimal  redesign  of 
existing  dynamic  structures  is  clearly  evident.  A 
good  example  is  aeroelastlc  tailoring  where 
stiffness,  mass  and  geometrical  characteristics  of 
an  aircraft  structure  are  optimally  distributed 
using  composite  materials,  for  the  improvement  of 
aeroelastlc  properties.  For  dynamic  structures  in 
general,  modifications  are  usually  focused  on  modal 
characteristics.  However,  since  modal 
characteristics  and  design  vectors  are  related  by 
eigenvalue  problems,  the  relationship  is  naturally 
complex  in  addition  to  being  numerically  formidable 
for  high  dimensioned  systems.  The  optimization 
problem  is  further  complicated  If  the  designer 
wishes  to  consider  a  large  number  of  constraints, 
mostly  inequalities,  on  the  minimum  and/or  maximum 
allowable  design  vector. 

In  this  paper,  the  optimal  redesign  problem  is 
treated  as  a  general  nonlinear  programming  problem. 
The  cost  function  and  constraint  equations  are 
linearized  about  a  nominal  design  to  obtain  a 
sequential  set  of  linear  optimization  problems.  The 
above  linearization  assumes  a  good  Initial  guess  of 
the  design  vector.  Unfortunately,  for  problems 
which  includes  eigenvalue  constraints  and  are 
structurally  complex,  good  initial  guess  of  the 
design  vector  may  be  impossible  to  obtain,  or  at 
best  depends  heavily  on  the  designer's  experience 
and  the  type  and  amount  of  redesign  desired.  The 
above  need  Tor  good  initial  guess  can  be  overcome; 
at  least  to  a  very  significant  degree,  by 
Introducing  the  continuation  method  of  handling 
constraint  equations.  This  amounts  to  replacing  the 
original  constraint  equations  by  a  sequential 
neighboring  set  of  constraint  equations.  Earlier 
applications  of  continuation  techniques  are  given  in 
[1  ,2]  where  closed-loop  eigenvalue  constraints  are 
considered  in  the  context  of  control  system  design. 
The  present  work  was  motivated  by  the  Introduction 
in  [3]  of  sequential  linear  programming  algorithm  to 
optimize  damper  locations  for  vibration  suppression. 
In  [3,<0,  the  sequential  linear  programming  approach 
combined  with  the  continuation  method  is  used  for 
the  design  of  control  systems.  The  development  in 
this  paper  closely  parallels  the  work  In  [ ti ]  and 
differs  primarily  In  the  optimlzatlon/de3ign 
objectives,  i.e.,  structural  redesign  Instead  of 
controls  design. 

To  demonstrate  the  proposed  approach  of  optimal 
structural  redesign,  a  Finite  Element  (FE)  model  of 
a  cantilever  beam  Is  considered.  The  thickness  at 
the  nodes  are  chosen  as  design  parameters  for  a 
constant  width,  rectangular  cross  section  bean.  The 
cost  function  to  be  minimized  Is  chosen  as  the  total 
mass  of  the  beam  while  selected  natural  frequencies 
are  gradually  driven  to  desired  values  subjected  to 
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geometrical  Inequality  constraints  on  the  maximum 
and  minimum  thickness  along  the  beam. 

TRANSFORMATION  OF  A  GENERAL  NONLINEAR 
PROGRAMMING  PROBLEM  TO  A  SEQUENTIAL  LIN¬ 
EAR  PROGRAM 


Applying  the  continuation  technique  of  Eq.  !3) 
to  the  right  hand  side  of  Eq.  (2),  we  obtain 
3  J 

maximize  J(p  )  *  — I  Ap  ♦  ...  (4) 

Ap  ip  **1-1 


Let  us  consider  the  general  nonlinear 
programming  problem, 


maximize 

P 

subject  to 


f(p)  K.  >}  f 


subject  to 
3f 

f(p.  .)♦  —  L  AP*. (1-Y.)f(pS)*Y  f° 
'  3p  Pl-1  1 


P.  *  Pi-, 


where  p  is  the  design  vector  and  f  Is  the 
constraint  objective.  By  linearizing  Eq.  (1 )  about 
a  nominal  design,  l.e.  at  1-th  step,  we  get, 


maximize 

AP 

subject  to 


J(P1.,)*— Ip  t 

3p  Pl-1 


f(p,_ .)  ♦  — L  Ap  ♦  ...  f° 

3P  p  i-1 


P1  "  pl-1  *  Ap 

The  formulation  In  Eq.  (2)  suffers  two  major 
problems;  a  feasible  solution  to  the  original 
problem  may  not  exist,  and  a  sufficiently  good 
Initial  guess  to  validate  linear  expansion  maybe 
impossible  to  obtain.  The  continuation  method 
resolves  the  above  problems  by  (1)  seeking  at  least 
a  feasible  solution  to  a  neighbor  of  the  original 
problem  If  Indeed  a  feasible  solution  to  the 
original  problem  does  not  exist  and,  (11) 
eliminating  the  need  for  a  good  Initial  guess  by 
starting  the  Iteration  with  a  neighboring  converged 
solution.  The  continuation  method  does  the  above  by 
replacing  the  original  constraint  objectives,  f°,  by 
a  sequential  neighboring  set  of  constraint 
objectives,  F(t),  where 

F( Y j )  -  (1  -  Yt)  f(p9)  ♦  Yj  r°  (3) 


where  p  Is  an  arbitrarily  chosen  starting  design 
vector  (more  appropriately,  p  Is  associated  with 
the  unmodified  existing  structure)  and  Y^  Is  a 
scalar  parameter  satisfying 


The  optimization  problem  of  Eq.  (4)  must  be 

transformed  to  a  non-negative  variable  so  that  the 
Simplex  algorithm  [5]  may  be  used  directly.  In 

addition,  the  output  of  a  Simplex  algorithm  may 
still  predict  large  corrections  for  some  elements  of 
Ap.  For  the  above  reasons,  we  Introduce  bounds  on 
the  maximum  allowable  local  parameter  change, 

-t  <.  AP  <.  t  (5) 

which  enhances  local  linearizations.  The 

transformation  to  non-negative  variables  and  the 

bounding  of  local  change  can  be  easily  done  by 
introducing  a  linear  transformation. 

y  -  Ap  ♦  e  (6) 

Using  Eq.  (5)  and  (6)  In  (4),  we  obtain  the  standard 
linear  program  at  step  "l”. 


maximize  — |  y  ♦  . . . 
y  3p  pl-1 


subject  to 


— L  y«.-.>)ti-Y1)f(ps)*Y.f0-f(p1.,)* 
3  p  pl-1 


y  £  2c 
y  >  0 


Ap  -  y  -  e 
P,  •  P.-i  *  Ap 


0  -  <  *2  •••  <  *n  *  ’ 

The  convex  combination  of  starting  constraint  and 
final  constraint  shown  In  Eq.  (3)  reveals  the 
following  facts: 

at  t  •  YQ  •  0,  F(0)  -  f(p3) 

at  Y  -  Yn  -  1 .  F(1  )  -  f(p°) 

l.e.  If  convergence  is  achieved  at  Y  -  1 ,  we  recover 
the  original  objective  constraint  condition.  It  Is 
also  obvious  but  nevertheless  important  to  note  that 
If  the  problem  has  a  feasible  solution  at  Y  -  1 ,  the 
problem  will  be  solved  In  N  steps. 


MINIMUM  WEIGHT  DESIGN  OF  A  BEAM  FINITE 
ELEMENT  MODEL 


To  demonstrate  the  approach  outlined  earlier,  we 
consider  here  the  problem  of  minimum  weight  design 
of  a  cantilever  beam  structure  modeled  by  finite 
elements.  A  nominal  uniform  beam  Is  assumed  given 
and  Its  corresponding  natural  frequencies  known. 
The  problem  Is  to  find  the  minimum  weight 
configuration  among  all  those  which  satisfies 
desired  natural  frequency  equality  constraints  and 
thickness  Inequality  constraints.  We  note  here  that 
the  above  constraints  may  arise  form  physical 
factors  such  as  material  property  limitations, 
desired  natural  frequency  locations  or  static 
buckling  design  limitations. 
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Consider  the  linear  free  vibration  equations 
M  x  (t)  ♦  K  x  (t)  -  0  (8) 

where 


subject  to 
Ey  <_  g  ♦  Ec 
y  <  2c 


31 


ne 

M  -  I  M 

j.l  J 


ne 

K  -  I  K. 
J-l  J 


(9s) 


—L  y  -  d-'»1)i(ps)*Y  x°-i(p1.1)  • 


represents  the  global  mass  and  stiffness  matrices. 
For  thin  beams  in  transverse  vibration,  the  j-th 
element  mass  and  stiffness  matrices  takes  the  form. 


1  D 

dp  V  1-1 

y  _>  0 

where 

ip  -  y  -  e 

Pi  ■  Pi-1  * 


(lib) 

(ne) 
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(9b) 
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u 

<J62 

d£2 

dC  (9c)  g  - 

P  -  P^, 

I 

-P  ♦  P^ 

* 

where  p  is  the  element  shape  function,  6.  and  n,  are 
the  mass  and  stiffness  distributions'3  over  Jj-th 
element  and  h  represents  the  element  length. 

To  form  the  design  vector,  we  let  the  thickness 
values  at  some  Judicious  choice  of  F£  nodes  be  the 
unknowns  and  linearly  interpolate  the  nodal  values 
whenever  we  need  Internal  values.  The  design  vector 
can  then  be  written  as 

P  •  <V  ••••  pnp)T 

where  P  represents  the  beam  thickness  at  1-th  node 
and  NP  Is  the  total  number  of  design  parameters. 

We  now  formulate  the  minimum  weight  eigenvalue 
placement  problem  as, 

minimize  J(p) 

P 


The  major  computational  effort  required  in  the 
problem  as  posed  in  Eq.  (11)  involves  solving  the 
generalized  symmetric  eigenvalue  problem 

XMv  -  Kv 

for  the  selected  NM  modes  and  the  solution  of  the 
linear  program  itself  using  the  Simplex  Method.  The 
eigenvalue  sensitivities  required  in  Eq.  (lid)  is 
derived  in  [6]  and  given  as. 
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-  v  (  —  - 

Ij  )  V 
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3Pj 

3PJ 

The  stiffness  and  mass  sensitivity  matrices  required 
above  can  be  computed  and  assembled  in  same  fashion 
as  stiffness  and  mass  matrices  themselves  In  Eq. 
(9) ,  and  is  given  by, 


where 

J(p)  -  pah[1/2,l,1,...,1,1,1/2]p 
subject  to 


(10a)  — 


AJ  ’  X  J 


Pl  <  P.  <  PU 

J  J  J 


1 . NH 

1 . NP 


(10b) 

(10c) 


where  X  is  the  objective  (desired)  eigenvalue  and 
p  ,  pu  are  the  lower  and  upper  bounds  on  nodal 
thicknesses. 

Following  earlier  derivations,  the  problem  posed 
in  Eq.  (10)  can  be  transformed  to  the  form  of  Eq. 
(7)  to  obtain  the  linear  program  at  step  1, 
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Ho  emphasize  here  that  In  the  above  problem,  the 
continuation  family  of  constraint  objectives,  as 
given  by  Eq.  (3),  was  applied  only  to  the  eigenvalue 
constraints  of  Eq.  (10b)  since  It  was  deemed  the 
only  coopromlsable  and  nontrivial  constraint. 


NUMERICAL  EXAMPLE 


Table  1  gives  the  data  for  the  particular 
cantilever  beam  FE  model  considered.  Two  cases  were 
computed  to  Illustrate  the  main  features  of  the 
design  approach,  specifically,  the  effect  of 
relaxing  constraints  and  the  effect  of  different 
starting  beam  configuration  on  the  final  converged 
configurations. 

The  effects  of  relaxing  stiffness  constraints 
are  given  In  Table  2.  The  starting  frequencies  and 
objective  natural  frequencies  are  given  In  Table  2a. 
The  total  mass  design  history  is  given  In  Table  2b. 
It  Is  noted  that  all  three  cases  reached  the  final 
objectives.  In  addition,  the  total  mass  decreases 
as  the  lower  bounds  on  thickness  decreases,  thus 
verifying  the  fact  that  more  design  flexibility 
leads  to  better  performance.  Table  2c  shows  the 
uniform  Initial  design  vector,  the  Imposed  lower  and 
upper  bounds  and  the  final  design  configurations, 
which  Is  highly  non-uniform.  Incidentally,  several 
elements  of  the  design  vector  have  reached  their 
lower  and  upper  limits  and  clearly  this  Is  useful 
Information  to  the  designer. 

The  effects  of  different  starting  beam 
configurations  are  give  In  Table  3.  Obviously, 
different  nominal  beams  corresponds  to  different 
Initial  natural  frequencies.  However,  all  three 
completely  different  starting  beams  converged  to 
essentially  the  same  total  weight  and  similar 
configuration  beam  as  seen  In  Tables  3b  and  3c.  The 
results  above  are  consistent  with  a  conclusion 
reached  In  [3]  that  the  continuation  procedure  with 
sequential  linear  optimization  more  often  converges 
and  yields  the  same  solutions  for  different  Initial 
conditions  than  does  conventional  nonlinear 
optimization  routines. 

In  the  previous  examples,  the  thickness 
constraints  and  the  objective  frequencies  were  not 
too  demanding.  The  above  relaxed  circumstance  were 
the  main  reasons  for  the  success  of  all  previous 
runs.  Additional  results,  not  shown  here,  Indicates 
that  Imposing  more  restrictive  bounds  on  the 
thickness  made  the  realization  of  frequency 
objectives  more  difficult  and  In  some  cases  its 
convergence  to  desired  conditions  impossible. 


CONCLUSIONS 


A  sequential  linear  programming  approach 
combined  with  a  continuation  method  of  handling 
constraints  has  been  derived  for  attacking  a  class 
of  nonlinear  programming  problems  and  In  particular, 
optimal  structural  redesign  problems.  The  numerical 
robustness  of  the  proposed  algorithm  has  been 
demonstrated  by  a  minimum  weight  redesign  problem 
with  eigenvalue  constraints  consisting  of  20  degrees 
of  freedom  and  38  constraints  with  11  design 
parameter. 


A  major  advantage  of  the  llncjr  optimization 
approach  presented  here  is  that  at  each  continuation 
step,  an  optimal  solution.  If  It  exists,  can  be 
found  very  efficiently  using  the  well  founded 
Simplex  method.  The  reason  Is  that  only  a  finite 
number  of  feasible  possibilities  exists  and  the  he 
Simplex  algorithm  efficiently  computes  the  optimal 
solution.  A  second  major  advantage  Is  the 
flexibility  to  handle  equality  and  Inequality 
constraints  on  both  the  design  variables  and 
functions  thereof. 

Finally,  and  perhaps  most  Importantly,  failures 
to  reach  the  final  (t  »  1)  solution  is  usually 
softened  by  convergence  to  an  Intermediate  (0  <  T  < 
1)  neighbor.  The  active  constraint  set  and  gradient 
Information  of  the  final  convergence  provides  a 
basis  for  Intelligent  revisions  of  the  problem 
statement. 
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Young's  Modulus 

21.5  X  106  psl 

mass  density 

.065  lb/in3 

beam  width 

1  in 

beam  length 

100  in 

number  of  elements  (uniform  length) 

10 

number  of  degrees  of  freedom 

20 

dimension  of  design  vector 

11 

number  of  equality  constraints  on  frequencies 

5 

number  of  lower  bound  Inequality  constraints  on  thickness 

11 

number  of  upper  bound  Inequality  constraints  on  thickness 

11 

number  of  inequality  constraints  on  maximum  allowable  local  parameter  change 

- 

11 

Table  1.  Data  for  graphite  epoxy  cantilever  beam  finite  element  model 


uS  (rad/sec) 

u°  (rad/sec) 

1.64 

2.82 

11.5 

10.9 

32.3 

30.8 

63.5 

67.0 

105.1 

104.8 

Table  2a.  Effect  of  relaxing 

conatraint s-star ting 
and  objective  frequencies. 

PS  -  (1 . 1)T 

pU  -  (2 . 2)T 


Y 

P*  *  -5 

P1  -  .3 

*  , 
p  -  .1 

0 

6.5 

6.5 

6.5 

.1 

6.1 

6.1 

6.1 

.2 

5.8 

5.8 

5.8 

.3 

5.7 

5.6 

5.6 

.4 

5.5 

5.2 

5.1 

.5 

5.6 

5.1 

4.5 

.6 

5.7 

5.1 

4.3 

.7 

5.8 

5.2 

4.4 

.8 

5.9 

5.3 

4.5 

.9 

6.0 

5.4 

4 .6 

2.0 

6.1 

5.6 

4.7 

Table  2b.  Effect  of  relaxing  constralnts-total 

mass  histories  for  different  stiffness 
constraints 


Table  3c.  Effect  of  different  starting  beam  configurations  -  starting  and  final  thickness 
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Abstract 

Several  methods  are  presented  for  place- 
ment/constralned  optimization  of  the  closed  loop 
eigenvalues  and  eigenvectors  of  linear  dynamical 
systems.  A  Unified  approach  Is  taken  to  (1) 
iterate  the  design  parameters  In  the  plant  being 
controlled,  (11)  ,the  location  of  sensor  and  actua¬ 
tors,  (111)  the  elements  of  a  direct  output  feed¬ 
back  gain  matrix,  and  (1v)  the  weight  matrices  In 
a  time-domain  LQG  performance  Index,  or  (v)  a 
combination  of  the  foregoing,  to  accomplish  a  con¬ 
strained,  simultaneous  optimization  of  the  sys¬ 
tem  s  closed  loop  eigenvalues,  eigenvectors,  and 
their  sensitivities.  A  low  dimensioned  discrete 
system  and  an  order  42  model  of  a  flexible 
structure  controlled  via  direct  output  feedback 
are  used  to  Illustrate  the  approach. 

1.  INTRODUCTION 

We  show  below  that  several  different  methods 
to  design  high  order  linear  feedback  control  laws 
can  be  unified  In  the  sense  that  a  single  approach 
can  be  taken  to  "optimally  tune"  these  methods 
vis-a-vis  the  placement  and  sensitivity  of  the 
closed-loop  eigenvalues  and  eigenvectors.  In 
Sections  2  and  3,  we  formulate  a  generalized 
closed  loop  eigenvalue  sensitivity  and  constrained 
optimization  approach.  In  Sections  4,  5,  6,  we 
show  how  three  different  approaches  to  design  of 
linear  feedback  controllers  lead  naturally  to  the 
problem  formulated  In  Sections  2  and  3.  In 
Section  7,  we  sunmarlze  numerical  solutions  of  two 
tximples,  finally  In  Section  8,  we  offer 
concluding  remarks. 

2.  CLOSED  LOOP  EIGENVALUE  PROBLEM 

We  are  concerned  with  linear  dynamical  sys¬ 
tems  In  the  generalized  state-space  form 

Zx  *  Ax  +  Bu  ♦  d  (1) 

with  1  inear  output 


and  linear  feedback  control 


H  Is  an  r  x  n  real  matrix  ! 

G  Is  an  m  x  r  real  control  gain  matrix 

and  the  closed  loop  system  matrix  Is 

I  -  A  ♦  BGH  (5) 

In  the  absence  of  disturbances  (d  *  0),  the  closed 
loop  performance  can  be  obtained  by  assuming  an 
exponential  solution  of  the  form  x  *  oexpD*); 
this  leads  to  the  generalized  eigenvalue  problems: 

right:  XjZa,  -  left:  x^e,  *  ATb,  (6) 

where  1  »  1,2,..., n  and  a,,  s,  are  the  right  and 
left  eigenvectors  corresponding  to  the  assumed 
distinct  eigenvalues  (x, . x  ).  The  eigen¬ 

values  and  eigenvectors  are  generally  complex.  We 
adopt  the  conventional  normalizations: 

[BlTZ(a)  -  |I)  ,  Ib]TA|o]  »  d 1 ag [ x , ]  (7) 

where  [a]  ■  [a,. ..a  I  and  [ b ]  *  (Bi-.-B.I  are  the 
right  and  left^'modSr  matrices.  1 

Let  us  now  address  the  situation  In  which 
ail,  or  at  least,  some  of  the  matrices  A,  Z,  B, 
H,  G,  and  therefore  A  •  A  +■  BGH,  are  functions  of 
a  q  x  1  system  design  vector  p.  The  elements  p 
Of  p  can  be,  for  example,  (1)  the  control  gains 
(elements  of  G),  (11)  an  Indirect  parameterization 
of  G  (e.g.,  the  weights  In  an  LQG  performance  mea¬ 
sure),  (111)  sensor/actuator  locations  (parame¬ 
terization  of  the  elements  of  B,  H),  or  ( i v )  plant 
model  parameters  (parameterization  of  the  elements 
of  A,  Z).  Since  A  «  £(p)  and  Z  -  'Z(p),  It  Is 
evident  that  x,  »  x^(p)  and  a,  *  o  (p):  except  for 
Isolated  events  (e.g.  blfunJatlort  points,  near¬ 
multiple  eigenvalues,  etc.),  we  can  consider 
x,(p)  and  o^(p)  to  be  continuous  and  differenti¬ 
able.  It  Is  therefore  reasonable  to  question 
whether  or  not  It  is  feasible  to  "tune”  p  to  solve 
a  constrained  optimization  of  x^p)  and  0i(p). 

The  first  and  second  order  sensitivities  of 
the  eigenvalues  are  derived  In  references  |1)  and 
[2],  these  are  as  follows: 


Thus  the  closed  loop  system  is  governed  by 

Zx  1  Ax  ♦  d  (4) 

where 

Z,  A  are  n  x  n  real  matrices 
x  Is  an  n  x  1  real  state  vector 

u  Is  an  m  x  1  real  control  vector 

d  Is  an  n  x  1  real  disturbance  vector 

y  Is  an  r  x  1  real  measured  output  vector 

B  Is  an  n  x  m  real  matrix 


-I,!!-  .  eT,^A_ 

»P»*Pm  1 ,pi,pm 


«  3ptJpm  ' 


81,P1m°1#1  jpJ 


■i8ipit' 
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"  ,B1P1mVkP1t*1  *  ,kP1w,1,1f>1t"k 

*  k.i  k.i 


where 

(9) 

and 

P1m 

.  jL  -  x,  i£- 

'  3pm  1  Jpm 

(10) 

3* 

»Pl 

»z  th  »2z 

*  >pt  *  *ptJpm  *  *p»*pm 

are  determined  from  direct  differentiation  of  the 
parametric  system  model  (which  obviously  must  be  a 
continuous  parameterization).  Analogous  develop¬ 
ment  lead  to  the  sensitivity  of  the  eigenvectors, 
and  the  above  simplify  considerably  for  the  most 
usual  case  that  Z  «  the  Identity  matrix.  Clearly 
the  above  formulations  suffer  singularities  near 
multiple  eigenvalues.  It  should  be  noted  that 
having  completed  the  solution  of  the  eigenvalue 
problem,  the  evaluation  of  the  eigenvalue  partial 
derivatives  represent  a  rather  modest  additional 
computational  expense. 


3.  EIGENVALUE/EIGENVECTOR  PLACEMENT/OPTIMIZATION 


One  popular  course  (In  attempting  to  place 
eigenvalues  for  multiple-input,  multiple-output 
systems,  MIMO)  Is  to  make  use  of  various  de¬ 
coupling  devices  to  map  the  problem  Into  a  family 
of  “pseudo-equivalent"  single-input,  single-output 
systems  and  thereby  render  the  eigenvalue 
placement  problem  trivial.  The  problem  with  this 
entire  class  of  approaches  lies  In  the  fact  that 
artificial,  physically  meaningless  constraints  are 
invariably  Introduced  which  often  lead  to  poor 
controller  designs  as  well  as  numerical 
difficulties;  see  Ref.  ( 3)  for  a  discussion  of 
these  Issues  and  the  recent  literature.  In  lieu 
of  attempting  to  eliminate  the  redundancy  In  an  ad 
hoc  fashion  (more  parameters  to  specify  than  the 
typical  number  of  eigenvalue  constraints),  we 
elect  to  follow  the  pattern  of  Refs.  1,  4-7  and 
exploit  the  redundancy  to  optimize  a  performance 
measure. 


Consider  the  following  optimization  pro- 

D’em:  We  wish 'to  minimize  a  criterion  function 

-  function  (Xj(p) . Xn(p),  •1(p),...,on(p).P) 

■  f(p)  (11) 

subject  to  a  vector  function  of  m.  equality  con- 
it’-alnts 

g  -  g(p)  -  0  (12) 

and  a  vector  function  of  m1e  Inequality  con¬ 
straints 

hj  s  h(p)  <  hu  (13) 

where  g0,  h  ,  hy  are  given  constant  vectors  de¬ 


fining  the  family  of  admissible  designs. 

first,  consider  only  the  equality  constraints 
of  Eq.  (12).  The  starting  Iterative  Ps»art  may 
result  in  large  violations  of  Eq.  12  an  dis  often 
too  far  from  a  zero  of  Eq.  (12)  to  permit  reliable 
convergence  using  a  generalized  Newton  algor¬ 
ithm.  Let  us  first  find  a  feasible  p  which  satis¬ 
fies  Eqs.  (12)  and  then  Eqs.  (13)1,  before  we  con¬ 
sider  the  Issue  of  optimization.  1  In  lieu  of  g0, 

we  Introduce  the  'portable*  objective  constraint 
vector  gp  such  that 

«P  -  T«0  MI  -  r)g(PsUrt)  (14) 

with  a  homotopy  or  continuation  parameter  y 
satisfying  0  s  y  s  1.  Replacing  g0  In  Eq.  (12)  by 
g.  of  Eq.  (14),  and  considering  p  •  p(y),  gives  a 
homotopy  family  of  equality  constraints 

H(p(r))  -  t90  ♦  (1  -  *)g(Pstlrt)  9(P(t))  '  0 

(15) 

Notice  the  boundary  conditions  on  Eq.  (15) 

at  y  -  0:  K(p(0))  -  9(PsUrt)  *  9<P(0))  *  0  (16) 

at  t  -  1:  H(p( 1 ) )  -  90  -  g(p(1))  -  0  (17) 

from  which  we  conclude  that  p(0)  *  p.t,rt  satis¬ 
fies  Eq.  (15)  for  y  *  0,  and  p(l).  If  it  can  be 
determined.  Is  a  feasible  solution  of  Eq.  (12), 
since  Eq.  (15)  -  Eq.  (12)  for  the  special  case 
Y  ■  1.  By  sweeping  y  Slowly  and  Iterating  on 
p(y)  to  satisfy  Eq.  (15),  we  can  always  Initiate 
each  Iteration  with  a  nearby  converged  solution 
and  thereby  very  nearly  guarantee  convergence 
(failure  will  occur  only  In  the  event  of  locally 
Singular  events  such  as  bifurcations,  turning 
points,  etc.,  In  which  other  remedial  action.  Ref. 
8,  can  be  taken). 

Preparing  to  develop  an  Iterative  process  for 
p(y)uwe  expand  Eq.  (15)  about  some  estimate  for 
P(yi)  ;  we  seek  a  correction  ap  such  that  Eq.  (15) 
will  be  satisfied  to  first  order  as  - 

H(p(y,)  ♦  ap)  -  H(p(Tl))  ♦  [f^lap  -  0  (18) 

where 

"(P^))  -  Y,g0  -  <l-T1)g(Pstart)  '  9(P(^»  ('9) 

I 9H/»p] 1 ]  -  -(99(p(t1))/»p(y1)I  (20) 

Since  the  dimension  m  of  H  Is  typically  much  less 
than  the  dimension  oT  p,  Eqs.  (18)  are  usually 
underdetermined.  Thus  a  criterion  must  be  intro¬ 
duced  to  select  a  particular  solution  for  ap.  We 
choose  to  minimize  the  correction  norm  apTap;  this 
serves  to  be  as  consistent  as  possible  with  the 
small  correction  assumption  Implicit  In  retaining 
linear  terms  and  also  Is  an  appropriate  design 
philosophy  (large  design  changes  are  qualitatively 
less  attractive  than  small  ones!).  Minimizing 


p(y,)  Is  estimated  from  a  converged  solution  for  a  neighboring  p(y j) - 


.',pTW/.p  subject  to  Eq.  (18)  gives 

JP  »  -W1l|gl1  |T  \p(yj)  (21) 

Nhere  w  Is  a  positive  definite  weight  matrix.  He 
use  Eq.  (21)  to  Iterate  for  p(tj)  In  •  Newton-like 
mannjr  by  recursively  computing 

P<M>new  '  ^Vold  +  Ap1 

until  convergence  Is  achieved,  for  each  of  a 
sequence  of  t’s 


{0  «  y0  <  Yj  <  ...  ym  »  1} 


Since  p(y,)  Is  Initiated  from  a  nearby  converged 
solution,  we  can  Insure  (by  choosing  Ay  suffi¬ 
ciently  small  In  y.  ■  y,  ,  ♦  ay)  that  the  elements 
of  H(p(y. ) )  are  arbitrarily  small;  failure  of  the 
correction  in  Eq.  (21)  will  occur  only  If  the 
Jacobian  [iH/ap|{]  becomes  locally  rank  deficient. 

In  Ref.  1,  we  show  how  to  generalize  the 
above  continuation  method  for  equality  constraints 
in  a  fashion  which  locally  Includes  the  active 
(violated)  subset  of  Inequality  constraints  and  a 
gradient  projection  decrement  of  the  performance 
criterion  J  of  Eq.  (11). 

The  resulting  algorithm  has  been  used  suc¬ 
cessfully  as  a  nonlinear  constrained  minimization 
a’gorlthm  and  Is  superior  to  conventional  con¬ 
strained  minimization  (nonlinear  programing) 
algorithms,  owing  to  the  use  of  continuation  to 
enhance  convergence.  As  Is  Shown  In  Ref.  8,  the 
domain  of  reliable  convergence  Is  vastly  larger 
than  conventional  Newton  and  quasi-Newton 
algorithms.  Alternatively,  we  show  In  Ref.  7  a 
method  to  make  the  local  corrections  via  a  linear 
or  quadratic  programming  algorithm  Imbedded  in  a 
similar  homotopy  family.  In  any  event,  the  use  of 
adaptive  continuation  and  homotopy  algorithms  to 
enhance  convergence  Is  a  very  Important  device;  It 
insures  that  failure  of  the  optimization  algorithm 
*111  not  be  due  to  a  failed  local  linearity 
assumption. 

4.  DIRECT  OUTPUT  FEEDBACK  CONTROL 

An  easy  to  state  MIMO  control  design  approach 
s  to  seek  the  m  x  r  elements  of  the  gain  matrix  G 
’n  the  closed  loop  system  matrix  of  Eq.  (5)  such 
that  the  eigenvalues  and  eigenvectors  of  A  ■  A  ♦ 
3Gh  solve  a  constrained  optimization  problem  of 
the  form  of  Eqs.  (11)-(13).  The  deslgp  vector  can 
easily  be  expanded  to  admit  sensor/actuator  loca¬ 
tions  and  plant  model  parameters.  Reference  1 
treats  this  problem  In  detail  and  Includes  numeri¬ 
cal  solutions  for  more  than  50  design  variables. 
:eferences  4-6  present  other,  lower  dimensioned 
applications  of  this  generic  type,  although  the 
numerical  methods  do  not  make  use  of  the  minimum 
correction  norm  and  continuation  methods  of  the 
foregoing  section.  However,  the  minimization.  In 
Jefs.  5  and  6  of  eigenvalue  placement  sensitivity, 
represent  significant  applications  of  this 
approach  to  design  of  robust  control  laws. 


5.  OPTIMAL  TUNING  OF  LQG  REGULATOR  WFIGHTS 

Consider  the  minimization  of  the  classl 
linear  quadratic  Gaussian  performance  Index 

J  ■  i  /  (xTQx  «•  uTRu)dt  (2 

1  0 

subject  to 

x  -  Ax  ♦  Bu  ,  (2 

The  necessary  conditions  for  minimizing  Eq.  (2 
subject  to  Eq.  (24)  are.  In  addition  to  Eq.  (74 
as  follows  (Ref.  9) 


i  »  -Qx  -  A 
u  -  -R'Vk 


where  i(t)  Is  an  n  x  1  co-state  (Lagrange  mult 
pller  vector).  If  we  seek  a  feedback  contro 
then  we  assume 


and  Eqs.  (24-(26)  give  the  matrix  Rlccatl  equatl 
for  K(t) 


R  +  KA  +  ATK  -  MBR'VjK  +  Q  *  0 


for  the  Infinite  upper  limit  of  Integration 
Eqs.  (24),  R  *  0,  and  Eq.  (28)  becomes  the  alg 
bralc  matrix  Rlccatl  equation  which,  if  the  syst 
(24)  Is  controllable,  can  be  solved  for  the  s 
metric  positive  matrix  K  using  Potter's  method 
Ref.  10.  Since  K  Is  constant,  using  (27)  In  (2 
gives  a  constant  gain  feedback  and  results  In 
closed  loop  system  of  the  form  of  Eq.  (4)  with  Z 
I  and  closed  loop  system  matrix 

A  ■  A  -  BR'Vk  (2 

The  dependence  of  J  upon  the  Initial  state  x( 
can  be  eliminated  by  considering  x(0)  to  be  un 
formly  distributed  on  the  unit  sphere.  Then 
lieu  of  minimizing  J,  we  minimize  the  expectati 
«(J);  this  leads  (Ref.  4)  to 


*(J)  -  \  trace  [Kt(x(0)xT(0)) ] 


If  we  assume  x(0)  to  be  taken  Independently  and 
be  scaled  appropriately,  t(x(0),  xT (0) )  is 
Identity  matrix,  so  we  can  replace  the  function 
minimization  of  Eq.  (23)  by  the  algebraic  requir 
ment  of  minimizing  the  trace  of  K.  Since  K 
dependent  upon  the  choice  of  weight  matrices  Q  aii 
R,  through  the  algebraic  solution  of  Eq.  (28),  « 
see  that  both  the  quadratic  performance  index  d 
Eq.  (30),  and  any  constraints  we  choose  to  impoj 
upon  the  closed  1m>(l  eigenvalues  and  elgenve 
tors  (of  A  ■  A  -  BR_1BtK)  depend  Implicitely  up 
Q  and  R.  The  question  of  choosing  Q  and  R  t 
achieve  eigenvalue  and  eigenvector  optimizatiq 
arises  quite  naturally,  and  the  problem  ha 
exactly  the  form  adressed  In  Section  3. 

Two  Issues  require  special  attention,  firs 
note  that  the  sensitivities  of  K  follow  from  dl 
ferentlatlon  of  Eq.  (28),  for  R  ■  0,  as  the  fol 


/.V.v'v:,'. 


•  V  /  -*  /  /  „•  ; 


lotting  Liapunov  equation 


S  ♦  V 


-  -  ♦  K  T-  IBR-1BT1K 


For  the  special  case  that  the  p  are  a  parame¬ 
terization  of  Q  and  R,  the  above  simplifies  to 

3p-I^*T|p_’*lp"‘  1(1881  »p~  r"1bTiic  (32) 

Secondly,  we  consider  the  fact  that  If  we  Iterate 
arbitrarily  upon  the  elements  of  Q  and  R,  they  may 
become  Indefinite  or  negative;  this  can  be  avoided 
by  using  Cholesky  decompositions 

0  •  [Q1/2I!Q1/2IT  .  R  -  (R1/2]IR1/21T  (33) 


0  ...  0 


q12  q22 


qln  q2n'*'  qr 


Thus  p  *  (qu  ...  qnn  ■;  rn  ...  r  1T,  and  using 
any  real  numbers,  for  q^  and  r^.will  guarantee  Q 
and  R  remain  positive  semi -definite. 

We  can  easily  augment  p  by  sensor/actuator 
locations,  etc.,  to  consider  more  general  elgen- 
*alue/e1genvector  optimization  problems. 

In  Ref.  11,  we  consider  the  more  general 
Quadratic  Index 


*jT}  [\  n]!> 

i  o  lui  In1  rj  lu) 


In  which  case  the  optimal  control  Is  given  by 

u  .  -R_1(BTK  ♦  NT)x  (37) 

and  K  satisfies  the  modified  Rlccatl  equation 

Q  ♦  *A  ♦  ATK  -  K|BR“1BTIK  -  0  (38) 

0*0-  NR*V  (39) 

A  *  A  -  BR'V  (40) 

ne  closed  loop  system  matrix  becomes 

A  *  A  -  BR‘lBTK  -  BR_1NT  (41) 

The  presence  of  the  cross-coupling  weight  N  *  0  In 


Cq.  (41)  obviously  affects  the  eigen-solution;  In 
Ref.  11,  we  establish  conclusively  that  Including 
K  allows  constructive  optimization  of  the  system 
eigenvalues  and  eigenvectors. 

6.  THE  GENERALIZED  PROBLEM  OF  MOERDER  AND  CALISE 

Reference  4  considers  the  problem  of  mini¬ 
mizing  c(J),  where 

J  *  i  J  (xTQx  ♦  uTRu)dt  ♦  f (G)  (42) 

0 

subject  to 

x  *  Ax  +  Bu  (43) 


Analogous  to  the  preceedlng  section,  c(J)  can  be 
written  algebraically  as 

t(J)  -  trace  [K]  +  f(G)  .  *{x(0)xT(0)}  *  I.  (46) 
where  K  satisfies  the  Rlccatl-llke  equation 

S(Q,R,A,B,G,K)  -  STK  +  KS  ♦  HTGTRGl||  =  0  (47) 

where  S  Is  given  by  Eq.  (5).  Reference  4  gives  an 
attractive  algorithm  for  Iterating  for  K,  G,  given 
Q,  R,  A,  B,  H.  Obviously,  the  design  vector  for 
this  approach  can  also  admit  parameterlzatlons  of 
Q,  R,  A,  B,  H,  and  algorithms  can  be  developed 
analogous  to  Sections  2-5  above  using  continuation 
methods,  as  necessary,  to  enhance  convergence. 

7.  NUMERICAL  EXAMPLES 

As  a  simple,  low-dimensioned  example,  con¬ 
sider  tuning  the  weights  for  an  LQG  controller  of 
the  type  In  Section  5,  with 


0  0  0  1  0  0 

0  0  0  0  1  0 

0  0  0  0  0  1 

-2  1  0  0  0  0 

1  -3  2  0  0  0 

0  1-10  0  0 


0  1/2 


P  M1  T 

T  "  LL 
|_N  RJ 


,p30”-p36  p8 


.1  ,  (Starting 


.1  ,  Estimate)  (50) 


We  adopted  the  following  performance  measure 
(eigenvalue  placement  penalty) 

1  3  2 

0(P)  ■  7  I  ♦;  h(*i)  (51) 

1  k-1  1  1 


•"**  *  ~t  ”  •  *  ,  «  , 


«,  ■  *,(P)  -  1  -  c,(p)  .  1  *  1.2.3  (52) 

h( • )  Is  the  heavyslde  unit  step  function 


damping  factor 


-Re(x1) 


m  Ml(p) 


(note,  the  eigenvalues  occur  In  three  complex 
conjugate  pairs  x.wRe(x.)±J  Im(x.)t  1  -  1,2,3). 
We  also  adopted  the  "three  equality  constraints 


91  “  "1o  *  “i 1  •  1.2,3  (54) 
and  of  course,  we  require  for  stability  that  c.  * 
0.  The  objective  values  of  «.  ,  were  simply  taken 
as  the  corresponding  zero  ga1™  values.  Observe 
that  the  performance  measure  of  Eq.  (51)  simply 
seeks  to  “herd"  the  eigenvalues  to  the  left, 
sub.ect  to  Eq.  (54),  which  attempts  to  maintain 
the  imaginary  components  constant  (this  Is  simply 
an  illustration).  In  Table  1,  the  progress  of  a 
continuation  process  Is  summarized.  The  homotopy 
was  designed  such  that  t  *  1  corresponds  to 
J(p(1))  *  0;  we  did  not  expect  to  be  able  to 
achieve  c.  *  1,  (l.e.,  drive  J  to  zero)  while 
holding  the  constant! 


Observe  in  Table  1  and  Figure  1  (the  corres¬ 
ponding  locus  of  the  eigenvalues  during  the 
continuation  Iterations)  that  as  t  was  swept  from 
zero  to  .90,  In  variable  steps,  the  eigenvalues 
marched  steadily  to  the  left  until  It  was 
impossible  to  maintain  the  Imaginary  parts 
constant.  Notice  the  dramatic  Increase  achieved 
in  damping;  not  all  weight  matrices  are  created 
equal!  The  Initial  diagonal  guess  (Eq.  (50))  on 
the  weight  matrix  square  root  converged  to 


.058  1.29* 


039  1.472 


-0.001  .101' 


00  00  -.001  .1012 

.043  -.202  -0.041  000  .174* 

-.177  .198  -0.026  000  -.013  .340Z 

L  (55)  ■ 

where  zero's  mean  the  converged  entries  were 
smaller  than  10'5.  Notice  that  the  lower  left 
comer  Is  non-zero,  thus  N  »  0  to  achieve  the 
solution  of  Table  1  and  Fig.  1.  Notice  In  Table  2 
xnd  Fig.  2  that  constraining  N  ■  0  for  this  case 
results  in  a  much  less  attractive  eigenvalue 
placement;  especially  for  the  third  "Inode.  Thus 
the  cross  coupling  terms  are  significant. 


As  a  second  example,  we  refer  tc  Ref.  1 
wherein  a  spatially  discretized  scheme  leads  to 
the  linear  second  order  system 


Mz  ♦  Cz  ♦  Kz 


to  model  a  flexible  structure  with  a  central  rigid 
hub  and  four  cantellvered  flexible  beams.  The 
dimensions  of  M,  C,  <  were  21  x  21;  the  first  10 
modes  were  considered  accurately  modeled.  A 


torque  actuator  acts  on  the  central  hub  and  two 
pairs  of  appendage  torquers  (nominally  acting  at 
the  mid  span  of  opposing  beams);  the  two  torquer 
pairs  are  constrained  to  Impart  Identical  torques 
to  the  structure,  so  effectively  u  Is  a  3  x  1 

vector  and  x  of  Eq.  (1)  Is  a  42  x  1  vector. 

Twenty  sensors  (ten  pairs  of  position  and  velocity 
sensors)  are  located  at  various  points  on  the 

structure.  The  above  can  be  cast  In  the  from  of 
the  problem  discussed  In  Section  2,  and  th 
algorithms  of  Section  3  Immediately  apply. 

Ref.  1,  we  solve  several  eigenvalue  optlmlzatlo 
problems  using  this  structure.  In  one,  w 

simultaneously  optimize  42  control  gains,  8  sensor 
locations,  and  2  plant  design  parameters  to  Impose 
12  eigenvalue  costralnts  and,  subject  to  these 
constraints,  minimize  the  sum  square  of  the 
control  gains.  The  continuation  process  of 
Section  3  converged  reliably,  even  starting  with 
zero  feedback  gains.  The  eigenvalue  trajectories 
during  the  continuation  Iterations  are  shown  In 
Fig.  3. 


8.  DISCUSSION 


* 


The  present  paper  establishes  Important  con¬ 
nections  between  various  approaches  for  designin 
constant  gain  linear  feedback  controllers,  an 
presents  a  unified  numerical  optimization  strateg 
for  simultaneous  or  sequential  tuning  of  control 
gains,  sensor/actuator  locations,  plant  desig 
parameters,  and  various  weight  matrices,  vis-a-vi 
constraints  and  optimality  criteria  specified  1 
terms  of  the  closed  loop  eigenvalues  an 
eigenvectors.  The  methods  presented  have  beenj 
found  rather  robust  with  respect  to  computational 
problems  one  usually  expects  for  high  dlmenslone 
nonlinear  optimization  applications. 


ill 


» 
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TABLE  1  Eigenvalues  for  LQG  Weight  Matrix 
Iterations  with  N  *  0 
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no  convergence  for  T  >  .9... 
TABLE  2  Eigenvalues  for  LQG  Weight  Matrix 
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Figure  3  DRAPER/RPL  Eigenvalue  Continuation  Locus 
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ABSTRACT 

A  mathematical  formulation  and  associated  algorithm  is  presented 
which  can  be  used  to  tune  the  weight  matrices  in  an  optimal  quadratic 
regulator  to  Impose  constraints  and  eigenspace  optimality  criteria 
upon  closed  loop  systems  eigenvalues.  The  algorithm  is  found  to  be 
efficiently  applicable  to  moderately  high  dimensioned  problems; 
reliable  convergence  has  been  routinely  demonstrated  with  over  one 
hundred  and  fifty  weight  matrix  elements  being  optimized  to  place 
eigenvalues  in  a  dynamical  system  of  order  fourteen.  These  results 
provide  a  basis  for  optimism  that  the  approach  is  applicable  to  a 
significant  family  of  problems. 

I.  Introduction 

The  design  of  practical  structural  control  systems  requires 
reliable  methods  to  determine  the  feedback  gain  matrix.  With  some 
approaches,  the  gain  matrix  is  directly  iterated  to  satisfy  design 
constraints,  while  with  some  others,  it  is  selected  by  solving  the 
Linear  Quadratic  Gaussian  (LQG)  regulator  problem. 

The  former  type  may  require  a  large  number  of  design  constraints, 
since  an  arbitrary  direct  feedback  gain  matrix  does  not  guarantee  the 
stability  of  closed  loop  system.  However,  the  latter  (LQG)  type  of 
design  may  encounter  other  difficulties,  as  examples,  the  resulting 
closed  loop  system  may  become  physically  meaningless  due  to  arbitrary 
weights  in  the  criterion,  or  numerical  difficulties  may  be  encountered 
in  solving  the  Riccatl  equation. 

In  recent  papers  by  Junkins,  et  al.,  (1,2),  structure  and  control 
design  techniques  to  satisfy  eigenvalue  constraints  have  been 


Introduced.  We  have  shown  that  the  direct  feedback  gain  parameters 
along  with  plant  parameters  can  be  iterated  via  a  nonlinear 
programming  method  based  upon  using  continuous  minimum  norm  correction 
strategy  and  homotopy  methods  to  enhance  consequence.  We  also 
Introduced  a  scheme  for  tuning  the  weight  matrices  in  a  LQG 
performance  criterion  to  achieve  eigenvalue  placement  constraints. 

This  paper  aims  at  developing  an  algorithm  for  optimally  adjusting  LQG 
weight  matrices  and  assessing  Its  effectiveness. 

In  Sections  2  and  3,  we  formulate  an  eigenvalue  sensitivity  and 
parameterization  scheme  considering  the  LQG  weight  matrices  as  design 
variables.  In  Section  4,  we  show  how  LQG  eigenvalue  placement 
optimization  can  be  formulated  as  a  nonlinear  programming  process.  In 
Section  5  and  6,  we  discuss  numerical  results  and  offer  concluding 
remarks. 

1 1 .  Closed-Loop  Eigenvalue  Problem 

Consider  the  linear  dynamical  system  in  the  state-space  form: 
x  =  Ax  +  Bu  (1) 

with  linear  feedback  control 

u  =  -Kx  (2) 

where  x  is  nxl  state  vector 

u  is  mxl  control  vector 

and  A,  B  and  K  are  plant,  controller  influence,  and  control  gain 
matrices  with  proper  dimensions.  We  assume,  for  initial  simplicity, 
that  the  full  state  is  measurable. 

Substituting  the  control  law  of  Eq.  (2)  in  Eq.  (1),  we  obtain  the 


closed-loop  system 


From  Eq.  (3),  we  have  for  the  present  case 


=  -  8 


The  evaluation  of  Eq.  (8)  depends  on  the  parameterization  scheme 
chosen  for  the  K  matrix,  which  will  be  discussed  in  the  following 
section. 


m 


III.  Regulator  Control  Problem 

In  this  section,  we  derive  optimality  conditions  for  three 
controllers:  direct  feedback,  LQG,  and  modified  LQG  types.  The 
closed-loop  stability  of  these  controllers  and  their  respective 
parameterization  schemes  are  discussed. 

The  derivatives  of  the  closed-loop  system  materix  (A)  required 
for  eigenvalue  sensitivity  calculation  are  derived.  Consider  the 
classical  LQG  problem: 


”  T  T 

minimize  J  =  J  (x  Qx  +  u  Ru)dt  (9) 

0 

subject  to  Eq.  (1)  where  the  weight  matrices  Q  and  R  are  assumed  to  be 
positive  definite. 

By  introducing  a  symmetric  positive  definite  matrix  Pss,  the 
optimal  solution  can  be  derived.  First,  we  rewrite  Eq.  (9)  as 

J  =  /  ( xTQx  +  uTRu  +  (xTP  x) ]dt  -  xTP  x  +  xTP  x  (10) 

g  dt  '  ss  '  =>  ss  =>  o  ss  o  v  ' 

where  x  and  x  are  the  state  values  evaluated  at  t  =  0  and  t  = 

o  tc 

respectively.  Assuming  that  our  closed-loop  system  is  asymptotically 
stable,  we  can  let  x  vanish. 
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Substituting  Eq.  (1)  Into  Eq.  (10),  we  then  obtain 
J  =  J  [xT(Q  +  P  A  +  ATP  )x  +  uTRu  +  2xTP  Bu]dt  +  xjp  x 

Q  ss  SS  SS  0  SS  0 

(11) 

In  this  equation,  we  substitute  the  direct  feedback  control  law  of  Eq. 
(2)  to  get 

J  »  J  xT(Q  +  Ps$A  +  ATPss  +  KTRK]xdt  +  xjp^  (12) 

where  A  Is  the  stability  matrix  In  Eq.  (3).  Given  Q,  R  and  K 
matrices,  the  stability  of  the  system  (3)  is  guaranteed  when  there 
exists  soltuion  P$s  to  the  Liapunov  equation  (71 

Q  *  PssA  ♦  ATPss  ♦  KTRK  -  0  '  (13) 

With  this,  we  can  formulate  a  parameter  optimization  problem  in  the 
form: 

minimize  J  =  xjpss*0  (14) 

K 

Subject  to  Eq.  (13).  If  we  consider  xQ  to  be  uniformly  distributed 
over  the  unit  sphere,  then  minimizing  the  expected  value  of  Eq.  (14) 

Is  equivalent  to  minimizing  the  trace  of  P$s. 

It  should  be  noted  that  stable  feedback  controls  requires  that 
all  the  eigenvalues  of  A  be  in  the  left  half  plane  of  the 
elgenspace.  Therefore,  if  we  attempt  to  impose  additional  constraints 
that  all  the  closed-loop  eigenvalues  have  negative  real  parts,  Eq. 

(13)  embodies  redundant  constraints. 

In  References  (1)  and  (2),  we  developed  an  algorithm  which  simply 
modifies  K  elements  without  the  constraints  (13)  to  enforce  all 
eigenvalues  being  in  the  left  half  plane.  Also,  we  can  consider 


generalized  performance  Index  of  Eq.  (14)  with  an  extra  penalty 
function  of  the  gains.  An  attractive  algorithm  of  Iterating  Ps$  and  K 
to  minimize  such  a  generalized  criterion  can  be  found  In  Ref.  [5J. 

Next,  minimizing  the  Integrand  of  Eq.  (11)  with  respect  to  K  and 
assuming  that  the  minimum  value  Is  zero,  we  obtain  the  classical  LQG 
optimal  control  law 

U  *  -  R_1BTPssx  (15) 


where  P$s  is  the  solution  of  the  algebraic  Riccati  equation 

g  ♦  PssA  ♦  ATP$s  -  P„«-Vp„  -  0. 

Thus,  the  minimum  performance  index  becomes  x 

J*  *  *IPssxo 


(16) 


(17) 


Note  that  If  the  open  loop  system  is  either  completely  controllable  or 
exponentially  stable,  then  the  solution  to  Eq.  (16)  exists  (7).  This 
condition  is  somewhat  less  restrictive  compared  to  the  case  of  direct 
feedback.  However,  we  have  to  solve  matrix  quadratic  equation  (16) 
Instead  of  linear  equation  (13)  for  Ps$.  Then,  as  long  as  either  the 
controllability  or  exponential  stability  of  open  loop  system  is 
maintained,  then  the  solution  can  be  found  by  using  Shur  method  [6]. 
Since  the  equation  (17)  is  a  Liapunov  stability  function,  the 
resulting  closed  loop  system  is  always  asymptotically  stable  (7). 

Now  we  consider  more  general  quadratic  index 


J  = 


(18) 


For  this  performance  index,  the  modified  optimal  control  law  and 
algebraic  Riccati  equation  are  given  by 


(19) 


U  =  -R-V  +  BTPss)x 

Q  ♦  PssA  ♦  -  P„«-Vp„  *  o 

where 

A  =  A  -  BR_1NT 


(20) 

(21) 
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Q  =  Q  -  NR_1NT  (22) 

Then,  the  corresponding  closed  loop  stability  matrix  becomes 

A  =  A  -  BR'W  +  BTPss)  (23) 

It  should  be  noted  that  cases  In  which  Q  is  not  positive  definite  or 
the  pair  (A,B)  is  not  controllable,  there  may  be  no  solution  to  Eq. 
(20).  But  if  a  positive  definite  solution  of  the  (Corresponding 
algebraic  Riccatl  equation  can  be  found,  this  solution  can  be  useful 
in  the  design  of  particular  control  system.  Allowing  the  "cross 
coupling  weight  matrix"  N  to  be  chosen  non-zero  is  shown  below  to 
permit  constructive  optimization  of  the  closed  loop  eigenvalues. 

Let's  assume  that  there  exists  a  positive  definite  solution  P$s 
although  the  theoretical  conditions  for  the  existence  of  P$s  (under 
these  generalized  circumstances)  have  not  been  clearly  defined.  Then, 
the  minimum  of  the  performance  index  can  be  found  and  becomes  Eq. 

(17).  Here  we  elect  to  parametrize  the  weight  matrix  in  Eq.  (18), 
which  must  be  symmetric  positive  definite.  Introducing  weight  matrix 
parameter  vectors  q,  r  and  n,  and  using  Cholesky  decomposition,  we 
rewrite  the  weight  matrix  in  the  form. 
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Therefore,  the  global  weight  parameter  vector  p  becomes 

P  *  Iqll*^21,‘** *qnn,rll*r21** * *rmm*nll*n21 . nmn,T  (26> 

For  the  calculation  of  closed  loop  eigenvalue  sensitivity  In  the 

previous  section,  we  need  to  differentiate  Q,  R  and  N  matrices  with 

respect  to  the  elements  of  p. 

The  partial  derivatives  of  Eq.  (24)  with  respect  to 

the  element  of  the  parameter  vector  p  leads  to 


f°T  N1  + 

.NT  rJ  3pa  3pi 


where  the  partial  derivative  of  L  can  be  obtained  by  direct 
differentiation  (as  simple  arrays  of  all  zeroes  except  a  unit  value  in 
the  element  corresponding  to  p^). 

Using  Eq.  (27),  we  write  the  partial  derivative  of  the  closed 
loop  stability  matrix 

f-  *  -  B  Hr  <mT  *  eTp5S>  -  BR'‘  <lr  *  8t  2s)  (28) 

*  *»  L  t 

It  can  be  shown  that  the  partial  derivatives  of  Pss  are  obtained  by 
solving  the  algebraic  Liapunov  equation 


SPss  A  +  AT  8Pss 


3P 


3P, 
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lA_  .  2/L  P  +  P  bR'Vp  (29) 

ss  3p4  sp^  ss  ss  ss  K  ‘ 


where  A,  Q  and  A  are  given  In  Eqs.  (21),  (22)  and  (23), 
respectively.  The  solution  of  Eq.  (29)  are  derived  In  Ref.  (1|  and  Is 
given  below. 

3P. 


Ipf  *  l.l'lT,!!.! 


(30) 


where 


'V.j 


iL((alT|l<HS1l|.||ij. 


(31) 


(a],  [ b ]  are  n  x  n  modal  matrices  whose  columns  ar^  a^,  b^, 
respectively  and  [RHS^]  Is  the  right-hand  side  of  Eq.  (29). 

Therefore,  using  Eq.  (28),  the  first  and  second  derivatives  of 
eigenvalues  In  Eq.  (6)  and  (7)  can  be  efficiently  calculated. 

For  the  classical  LQ6  case  of  N  *  0,  the  Q  and  R  matrices  may  be 
parameterized  separately  by  using  the  same  scheme  presented  here  with 
q  and  r  vectors. 

IV.  Optimization  Approach  for  Eigenvalue  Placement 

The  problem  of  eigenvalue  placement  can  be  formulated  by  an 
optimization  approach  In  which  we  seek  to  impose  specified  eigenspace 
constraints  and  minimize  an  eigenvalue  placement,  sensitivity  or  other 
closed  loop  eigenvalue/eigenvector  performance  criterion. 

Utilizing  the  parameterization  scheme  described  in  the  previous 
section,  we  consider  the  following  sets  of  equality  and  inequality 
constraints: 


*oi 


g<(p)  =  o 


1  —  1,2,... m. 


(32) 
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foj  •  f j(p)  >  0  j»1.2....m1e  (33) 

where  g0^  and  are  objective  constraint  function  values  or 
constraint  boundaries  corresponding  to  the  1th  closed  loop  eigenvalue, 
and  g^(p)  and  f^(p)  are  current  values  of  these  constraint 
functions.  The  feasible  solution  to  Eqs.  (32),  (33)  can  be  obtained 
by  either  considering  the  locally  active  Inequality  constraints  as 
equality  type,  or  by  minimizing  a  bounded  quadratic  penalty  function 
subject  to  the  equality  constraints.  We  choose  the  latter  approach, 
which  leads  to  an  optimization  (nonlinear  programming)  problem  of  the 
form: 


1  2 

minimize  4  £  K.^h(-4>.) 

Li  11  I 


(34) 


subject  to  Eq.  (32)  where 
♦l  ■  fo1  - 

K1  >  0  is  a  weighting  factor  and  h(-)  is  the  heaviside  unit  step 
function.  This  nonlinear  programming  problem  can  be  solved  by  using 
continuous  minimum  norm  correction  algorithm  of  Ref.  (3).  The 
essential  feature  of  this  algorithm  is  to  solve  a  set  of  nonlinear 
equations  for  variable  continuation  steps.  Its  first  step  is  to 
generate  a  homotopy  family  of  problems  by  introducing  "portable" 
objective  function  values  defined  by  the  linear  map 

9P( y )  =  y90  +  (1  -  Y)g(PsUrt)  (35) 

with  a  homotopy  or  continuation  parameter  y  satisfying  0  <  y  <  1.  If 
we  replace  the  original  objective  g0  in  Eq.  (32)  by  the  portable 
objective  gp(y),  we  obtain  a  homotopy  family  of  constraint  functions: 
H(p(r))  =  y90  +  (1  -  r)g(Pstart)  -  g(P(v))  =  0  (36) 


Thus,  at  y  *  0,  the  above  has  the  solution  Pstart*  an<*  y  -  1.  it 
becomes  the  original  equations,  so  p(l)  Is  the  desired  solution. 

Since  Eq.  (36)  is  usually  underdetermined,  a  unique  correction 
vector  &p  can  be  obtained  by  a  minimum  norm  solution  of  the  truncated 
Taylor  series  expansion  of  Eq.  (36):  that  is,  we  minimize  apTap 
subject  to  the  equality  constraint 

H(p(y)  +  aP)  =  H(p(  y)  )  +  t-f^lAp  =  0  (37) 

Then,  the  minimum  norm  correction  vector,  ap  takes  the  standard  form 

*p-  0Ti'lp|1lp|Trl4H  (38) 

where 

AH  =  -H(p( y) ) 
and 


[iHi  =  _  tag(p) | 

ap 1  1  ap  1 


Using  Eq.  (38)  and  starting  with  a  neighboring  solution  at  step  y  = 
y i _i ,  we  refine  p(y^)  recursively  by  computing 

P<Vnew  =  P(’l>old  *  (39) 

until  local  convergence  for  each  y^J  final  convergence  is  achieved  by 
incrementing  y^  after  each  local  convergence  until  y  approaches  1  (or, 
if  y  *  1  cannot  be  achieved,  continue  the  process  until  y  is  as  large 
as  possible,  i.e.,  we  accept  the  final  local  convergence).  Along  with 
the  Iterations  to  satisfy  the  equality  constraints,  the  performance 
index  can  be  reduced  by  introducing  the  performance  index  as  an  extra 
equality  constraint  with  zero  objective  value.  Incidentally,  this 
procedure  has  been  shown  theoretically  equivalent  to  the  gradient 


projection  algorithm  for  constrained  optimization  problems  [9|,  If  an 
equivalent  homotopy  Imbedding  Is  Introduced  to  control  step  size  In 
the  classical  gradient  projection  method. 


V.  Computational  Study 

As  a  test  example  to  demonstrate  optimal  tuning  of  the  weights 
for  an  LQG  controller,  we  selected  the  DRAPER/RPL  model  [1].  It 
consists  of  a  central  rigid  hub  and  four  arms  (identical  cantilever 
beams)  with  tip  masses;  In  addition.  It  is  assumed  that  torque 
actuators  are  located  at  the  center  of  the  hub  and  at  the  middle  span 
of  each  appendage.  The  model  is  described  by  Eq.  (1)  in  which  the 
order  of  the  system  Is  n  =  14  and  and  the  number  of  actuators  ism* 

3.  The  number  of  controller  inputs  is  three  (instead  of  five)  since 
the  actuators  or  opposing  appendages  are  constrained  to  apply 
Identical  control  torques.  The  detailed  configuration  and  the  nominal 
structural  parameter  values  may  be  found  in  (1). 

For  this  example,  we  adopted  three  equality  and  five  inequality 
constraints  on  the  dominant  closed  loop  eigenvalues,  i.e. 


uoi  -  wi<P)  =  0  1  =  1.2.3  (40) 

1  -  ^(p)  >0  j  =  1,2,. ...5  (41) 

where  and  are  undamped  frequency  and  damping  factor 
corresponding  to  the  i^  eigenvalue,  respectively.  We  implicitly 
require  that  all  remaining  eigenvalues  remain  in  the  stable  left-half 
plane.  The  objective  frequency  values  were  taken  as 


uq1  =  «53  ,  uj_ 


=  4.38 


7.91. 
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The  equality  constraints  of  Eq.  (40)  were  employed,  for  simplicity  and 

to  avoid  possible  bifurcations  leading  to  unnecessary  numerical 


complications,  so  that  the  desired  eigenvalue  locus  will  remain 
horizontal.  The  objective  damping  factors  were  taken  as  the  critical 
damping  condition  (cQ^  =  1);  these  were  chosen  to  move  the  eigenvalues 
as  far  as  possible  to  the  left-hand  side  (these  critical  damping 
factors  obviously  cannot  be  achieved  while  holding  the  frequencies 
constant!).  However,  progress  toward  this  objective  will  tend  to 


Increase  the  damping  factors  as  much  as  possible  and  will  be 
constructive  vis-a-vis  eigenvalue  placement;  we  Implicitly  expect 
final  local  convergence  for  some  y  <  1.  The  quadratic  penalty  terms 
for  each  mode  In  Eq.  (34)  were  equally  weighted,  i.e.,  =  1,  i  = 

1,2,..., 5.  The  total  set  of  8  constraints  of  Eqs.,(40)  and  (41)  were 
Imposed,  by  iterating  153  parameters  representing  Q,  R  and  N  matrices 


with  the  starting  diagonal  matrix 

L  =  diag  {1Z,12,12,12,12,12,12,.32,.32,.32,.32,.32,.32,.32,12,12,12} 


At  each  continuation  step  starting  with  y„  s  0,  y1  s  .1  and  ir  *  .1, 
the  convergence  to  the  intermediate  solutions  was  reliable  for  small 
values  of  y.  but  it  became  necessary  to  make  ay  adaptive  for 
larger  y  values.  If  the  local  iterations  at  y^+j  did  not  converge, 
then  the  step  size  Ay  was  reduced  by  half  until  it  became  less  than 
.005  without  achieving  local  convergence. 

In  Tables  1  and  2,  the  results  of  two  cases  1)  with  N  *  0  and  2) 
with  N  =  0  are  reported  for  variable  continuation  steps.  For  both 
cases,  the  damping  factors  for  all  modes  including  the  unconstrained 
modes  increased  dramatically  and  similar  eigenvalue  locus  were 
obtained  as  can  be  seen  in  Fig.  1. 

At  each  iteration,  the  norm  of  both  control  gain  vector  (K 
matrix)  and  eigenvalue  sensitivities  with  respect  to  gain  elements 


were  calculated.  Figure  2  and  Fig.  3  show  that  the  gain  and 
sensitivity  norms  for  both  cases  Increased  exponentially.  However, 
the  gain  norm  for  Case  1  was  much  less  than  that  for  Case  2,  and 
similar  result  for  eigenvalue  sensitivity  norm  was  obtained.  The 
number  of  gradient  calculations  for  the  optimization  were  counted  133 
and  99  for  Cases  1  and  2,  respectively.  Most  of  computation  time  was 
spent  In  solving  the  algebraic  Rlccatl  equation  by  Shur  method. 

Instead  of  Shur  method,  a  Newton  type  of  algorithm  [8)  can  be 
used  to  solve  the  Rlccatl  equation  and  to  (possibly)  save  computation 
time  (after  the  first  nonzero  y  continuation  step  is  made  by  the  Shur 
method).  The  convergence  of  a  Newton  algorithm  wip  likely  be 
quadratic  since  the  iteration  can  be  initiated  with  "arbitrarily 
close"  estimates  of  the  solution  vector  at  each  continuation  step. 
However,  the  modest  increase  in  storage  requirments  should  be 
considered. 

For  both  cases,  the  starting  diagonal  matrices  resulted  in 
convergence  to  fully  populated  weight  matrices  at  the  end  of  the 
continuation  Iterations.  Obviously  the  eigensolution  and  gain  matrix 
are  constructively  affected  by  properly  chosen  off-diagonal  elements 
(N)  of  the  weight  matrix. 

VI .  Conclusions 

This  study  presents  numerical  results  obtained  by  an  algorithm 
for  sequential  tuning  of  LQG  weight  matrices.  It  is  obvious  that 
physical  performance  measures  usually  does  not  dictate  unique  choices 
for  the  LQG  weight  matrices.  We  have  shown  that  the  fully  populated 
weight  matrices,  especially  Including  the  cross-coupling  terms  in  LQG 


performance  criterion,  affect  both  the  control  gain  and  closed  loop 
elgensolutlon.  Since  these  matrices  are  often  selected  as  simple 
diagonal  or  block  diagonal  matrices,  the  desirability  of  systematic 
methods  to  optimize  these  weights  is  evident.  The  results  of  the 
present  study  are  a  basis  for  optimism  that  for  systematic  eigenvalue 
placement  can  be  achieved;  the  present  method  has  been  found  reliable 
and  numerically  stable  for  the  test  examples  considered  In  this  study. 
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